Continuous Phase Modulation The wikipedia article on continuous-phase modulation states the following:

  For instance, with QPSK the carrier instantaneously jumps from a sine to a cosine (i.e. a 90 degree phase shift) whenever one of the two message bits of the current symbol differs from the two message bits of the previous symbol. This discontinuity requires a relatively large percentage of the power to occur outside of the intended band (e.g., high fractional out-of-band power), leading to poor spectral efficiency.

I want to know how this is true since in M-ary PSK the amplitude is constant so isn't the power constant too?

What is important is the envelope of the modulated signal. For M-PSK, due to the phase discontinuities, the envelope droops every time there is a phase transition. The envelope even goes to zero if there is a phase transition of $180$ degrees. So M-PSK does not have a constant envelope.

Changes in the envelope are largely reduced by using offset QPSK (OQPSK), and they are completely eliminated by using continuous-phase modulation schemes, such as minimum-shift keying (MSK), which is a special case of OQPSK.

Take a look at these related questions and its answers: Q1, Q2, Q3.


Energy of an audio signal When energy is referred to as the square summable sequence, is that energy a representation of the power of the signal? How is the energy of frequency determined?

E.g. what if an audio signal has high energy on the higher frequencies and another has high energy on the lower, bass, frequencies, how is this determined?

What exactly is the power from the square summable series representing then and is it useful for audio?

No, energy is not the same as power.  Power is energy per unit of time.  An analogy- energy is like distance, while power is like speed.

The easiest way to determine energy in the frequency domain is to take the signal and do an FFT.  The Fourier Transform preserves energy, i.e. the transform has the same exact energy as the signal does in the time domain.  You can then find the energy of each frequency bin by multiplying its value with its complex conjugate.

Signal power represents how quickly it delivers energy.  In the context of audio it correlates with how loud the sound is, though there is not a 1-to-1 correspondence.


Adding Dirac Delta with Unit Step What does the graph of the function $x(t)=-\delta(t)+u(t)$ look like?
$$\delta(t)\ldots\text{ Dirac delta impulse}\\
u(t)\ldots\text{unit     step function}$$

Will the impulse at the origin start from 0 downwards or from 1?
I'm having trouble because the unit step is discontinuous at t=0, so I don't know how it will shift the impulse.
I am talking about continuous time.

Note that $\delta(t)$ is no function but a distribution. This means that it doesn't have function values but it is only defined by its integral. Consequently, the (generalized) function $x(t)=-\delta(t)+u(t)$ has no value at $t=0$ (the only useful value would be $-\infty$ but that's kind of hard to sketch). For sketching the function graph, you can define function values for $t\neq 0$ in the usual way. At $t=0$ you simply sketch a negative $\delta$-impulse, for $t<0$ the function is zero, and for $t>0$ you assign it the value $1$.


Sampling frequency lower than half of sampled signal frequency? In GNURADIO when RTL-SDR Source block is placed there is a parameter Sample Rate (sps) and Channel frequency.

I can put there for example 32kHz into Sample rate and 107.5MHz into Channel frequency to tune some FM station in my location. FFT block shows a signal of that station in a correct frequency as shown on the picture:

But according to Nyquist–Shannon sampling theorem, Sample rate should be at least 2 times higher than signal frequency, so in this situation I should set at least 215MHz into Sample rate parameter.

Why can I set much lower sampling frequency of a signal and it still looks correctly?

The channel frequency controls the local oscillator on your SDR which is the frequency about which it covers.  So you are receiving signals from 16KHz below 107.5MHz to 16KHz above.

The concept you want to look up is called heterodyning.  When you modulate (multiply) a signal by a sine wave you end up with two copies shifted in frequency space by the frequency of the sine wave.  So in your case, you've told the local oscillator to run at 107.5MHz which shifts the 107.5MHz carrier from the radio tower to 0Hz and to 215MHz though the high frequency version gets cut off by the low pass filter in your radio.

If you are interested in the math then if we want to shift a complex signal $X(t)$ by a frequency of $f$ we multiply by $e^{2\pi f i t}$.  $i\sin(2\pi f t) = e^{2\pi f i t } - e^{-2\pi f i t}$ where $i = \sqrt{-1}$.  So multiplying by $i\sin(2\pi f t)$ gives $X(t)i\sin(2\pi f t) = e^{2\pi f i t}X(t) - e^{-2\pi f i t}X(t)$ which is where the two copies come from.  (The multiplication by $i$ in, e.g., $i\sin(2\pi f t)$ ultimately just adds a $90^\circ$ phase shift.)


Resampling and removing high frequency noise? I am currently working on a simple sampler that will allow me to load in a wav file and use my MIDI keyboard to play the loaded wav sample at the frequency according to the note played.

Now I need this pitch shifting to be done in real time so I opted to use resampleing to achieve different pitches. I don't care that the sample speeds up or slows down when doing so.

What I am doing is say for example sake I have a wav file of a sine wave at 261.63hz (middle C) and this file has 1000 samples (example). Now I want to resample this sine wave so it plays at 146.83hz. I am dividing the desired frequency of 146.83hz by 440hz which gives me 0.3337. I then divide up the 1000 samples by 0.3337 to get the desired frequency. I will admit I am fairly new to all this but I think I am doing this correctly? This now leaves me with a smaller array of samples that when played gives me a rougher 146.83hz sine wave.

Now the issue is I am getting high frequency noise in the 2k to 20k range which varies on each note played. I think I am suppose to use interpolation to remedy this but I have not figured out how to implement it properly in my code. Do I resample then interpolate the resampled data OR do I resample and add in the interpolated samples between each sample?

I tried linear, cubic and hermite interpolation but neither made a difference in removing the high frequency noise. Though I don't believe I am using them properly. I got the interpolation the code from http://www.musicdsp.org so I'm  pretty sure it's correct.

Any help or suggestions would be great.

Sample playback

The basic idea of sample playback in musical applications is to keep track of each voice's playback position, to form an output sample by reading the source sample data at the playback position, to add a possibly time-varying playback step to the playback position, and to repeat this in a program loop until we have accumulated enough output samples for that voice. For multiple voices, their output samples are summed in a mixing buffer. If the output sampling frequency differs from the source sampling frequency, or if we want to play a different note than the source note, then playback step may be not equal to 1 (one source sample per one output sample), using your example scenario:

                source sampling frequency * desired note frequency
playback step = --------------------------------------------------
                output sampling frequency * source note frequency

                44100 Hz * 146.83 Hz
              = -------------------- = 0.56121239919
                44100 Hz * 261.63 Hz

If the playback step is not an integer, then the playback position will sometimes also be not an integer and we need to interpolate between the samples. It is also possible to do vibrato and tone portamento by varying the desired note frequency with time and consequently also varying the playback step, which will inevitably lead to a mostly non-integer playback position and the need for arbitrary sample rate conversion.

Interpolation, spectral images and aliasing

In this context, sampling is best understood (see Andreas Franck (2012) Efficient Algorithms for Arbitrary Sample Rate Conversion with Application to Wave Field Synthesis, section 3.4.1 Sample Rate Conversion as an Analog Resampling Process) as multiplication of a band-limited continuous-time audio signal by a Dirac comb and storing the integral over each scaled Dirac pulse in the resulting impulse train as the discrete-time sample data. The Fourier transform of the impulse train is periodic with a period of $2\pi$ in angular frequency. The chosen interpolation method is characterized by its continuous-time impulse response, and interpolation can be seen as convolution (filtering) of the impulse train by the continuous-time impulse response. The impulse response of a reasonable interpolation method is approximately low-pass, meaning that the convolution attenuates the spectral images (the duplicate spectra centered at multiples of $2\pi$). Then sampling of the approximately low-pass signal, at the playback positions encountered, will hopefully not lead to significant aliasing of the spectral images into the audio band.

If the bandwidth of the audio signal is half its sampling frequency, then a playback step > 1 will lead to aliasing of not only the spectral images, but also of some of the base-band audio frequencies. It helps to have an intermediate output sampling frequency where there is headroom for aliasing into a band above the audio frequencies that can be filtered away at an additional sample rate conversion to the final output sampling frequency. It is also possible to have multiple versions of the sample data that have been pre-filtered to reduce their bandwidth, to avoid aliasing when playback step > 1, similar to MIP mapping in computer graphics, see Laurent de Soras (2005) The Quest For The Perfect Resampler.

A time-varying playback step may further increase the continuous-time bandwidth, but usually the variation is slow enough that the playback step can be considered locally constant when it comes to bandwidth considerations.

Looping

Traditionally samplers have a number of ways of playing back sample data:

no loop (one-shot),
forward loop, with integer loop start and integer loop length (in sampling periods equal to 1 / sampling frequency, or the time difference between successive samples), and
ping-pong loop (forward-backward-...), with integer loop start and integer loop length (in one direction).

The forward loop requires a change in the program loop. If the playback position goes to the loop end (the first sample not included in the loop) or beyond it, then the loop length is subtracted from it. This works assuming that the playback step is smaller than the loop length. Ping-pong loop requires a bit more program logic, or the ping-pong loop can be unrolled into a forward loop. Upon note release, it is also possible to let the playback position escape the loop and to play the source sample data until the end.

Considering forward loops, when playing back the source sample data at its sampling frequency, what the loop length being an integer number of sampling periods means is that the period of the output waveform will be quantized to an integer multiple of the source sampling period. So it typically won't be possible to produce exactly the correct frequency, unless the source note frequency is a submultiple of the source sampling frequency. With your example, a loop length of 1000 sampling periods gives 261.63 Hz * 1000 / 44100 Hz = 5.932653 periods of the sine wave, so the last period will be cut short and give you a wide-band buzzing noise no matter what interpolation method you use.

Because of this inflexibility, we can bring into question whether the loop length should be an integer number of sampling periods. The benefit of a floating point (or a fixed point) loop length would be that we get a more correct note frequency, but the danger is that we introduce a discontinuity in the interpolated signal. With an integer loop length this discontinuity can be avoided by doing a short cross-fade of the loop start and end so that they have a few identical samples around them, as many as needed by the interpolation method. This way the interpolated curve will be exactly the same at the beginning and at the end of the loop. An alternative to cross-fading is to duplicate a few samples from the beginning of the loop at the loop end and in the following samples and to move the the loop slightly forward in time so that the interpolation method sees the same samples at the start and at the end of the loop. In case of a floating point loop start and loop length, they can be adjusted until no discontinuity can be heard, although this ties high loop quality to the present choice of the interpolation method.


How is Wi-Fi transmitter sends complex number to receiver? I am recently studying physical layer the Wi-Fi (802.11), after doing the FFT, the signal bits are transformed from '0' and '1' to complex numbers, but the complex number is not exist in real life, so my question is:

How is transmitter send those complex numbers?
What is the transmitter actually send?

How is transmitter send those complex numbers?

It doesn't.  As you noted, that's impossible.

What is the transmitter actually send?

The short answer is that the transmitter takes a baseband signal expressed as inphase/quadrature (I/Q), and modulates the inphase part with $\cos \omega t$ and the quadrature part as $\sin \omega t$.  On the receiving end, the inphase part is demodulated by $\cos \omega t$, and the quadrature part is demodulated by $-\sin \omega t$.
This acts as if the I/Q signal were complex with the inphase part being real and the quadrature part being imaginary, and that the transmitter modulated it with $e^{i \omega t}$, and then only transmitted the real part.
It turns out that you can work out all the inphase/quadrature math, very laboriously using only real numbers and trig identities -- or, you can note that if you pretend that your I/Q signal is complex, then you can do all the arithmetic with complex numbers and everything behaves exactly the same.


LPF in the stage of IQ demodulator is it a analgor filter or digital filter? I recently designed the LPF of the IQ demodulator using the Butterworth LPF refering to https://dspillustrations.com/pages/posts/misc/baseband-up-and-downconversion-and-iq-modulation.html. But I have a question. If you look at the block in the picture below, you can see the baseband block after passing the A/D Conversion. I think I applied IIR digital filter, but the composition of the picture below shows that LPF is analog domain. So I'm suddenly very curious. Is this LPF analog filter or a digital filter?

I designed Lowpass filter corresponding to LP block using Butterworth IIR filter design method in the IQ demodulator block diagram above. I first calculated analog butterworth filter coefficient, and then I calculated IIR filter coefficient b,a of z-domain using the bilinear transform. And with this, we calculated the waveform through the lpf using the filter(b,a,x) function of the matlab. So, we got the spectrum for Inphase and Quadrature phase as below. But in this process, I have a question. Obviously, I designed a digital filter of z-domain, and I wonder why digital filter exists in analog domain because it is before ADC in block diagram.

The A/D can be placed as a single real A/D before the multipliers, OR as shown in the diagram as two A/Ds one after each multiplier to sample the I and Q channels. In either case, an analog filter is required before any A/D conversion as an anti-alias filter. This can be a bandpass filter or a low-pass filter, depending on which image in the analog domain's frequency spectrum is desired as the signal of interest. Even if there is no actual additional signal in the image bands, with no filter (and gain in the receiver) the noise figure will be significantly degraded due to the accumulated noise floor from each of the images. An analog filter prior to any sampling is important.
So in this case given the placement of the A/D converter, the filters as shown would be the analog anti-alias filter that I describe.
If the reason for image rejection is confusing, please refer to this posts that details how aliasing occurs in the digital sampling process: Higher order harmonics during sampling
So in summary, for the implementation of an IQ demodulator, there could be two A/D converters if I and Q are to be sampled separately: one for the I channel and one for the Q channel resulting in a complex representation of the baseband signal. Alternatively there can be a single A/D converter prior to the multipliers and then the multiplier and filter structures as shown can be implemented all digitally (referred to as a Digital Down-Converter or DDC). In this case, a single analog anti-alias filter (band-pass or low pass as desired for the implementation) is still needed in front of the A/D converter, in addition to the digital low pass filters after the multipliers to select the low frequency component after the multipliers.
The block diagram as shown in the question does not make that clear that $v(t)$, $y(t)$ and $\hat d[k]$ are complex signals (although clearer if you read through the linked reference).


generating white gaussian noise in matlab using two different functions I want to know the difference between the two Gaussian noises generated below? Which one is white and how can i make the other one white?

y=wgn(1,10000,0)

and

y=randn(1,10000);

A white noise sequence is one for which each (random) element is uncorrelated from every other element:
$$
E[y[n]y[m]] = \left \{ \begin{array}{ll}
0 & \mbox{for } n\not=m\\
\sigma_y^2 & \mbox{for } n = m
\end{array}
\right . \\= \sigma_y^2 \delta[n-m]
$$
where $\sigma^2_y$ is the variance of $y$.

Note that I am assuming (because of the whiteness) that the signal is zero mean.

The whiteness of a signal says nothing about the distribution of its values.  To know something about that, Gaussianity or some other distribution needs to be invoked.

The functions wgn and randn both produce white, Gaussian noise sequences.

Calling the function rand would produce a white, uniformly distributed noise sequence.


Can use of Fourier transform be minimized completely with the help of Laplace and Z transform? Fourier transform has different types like continuous Fourier transform (CFT), Discrete time Fourier transform (DTFT) and Discrete Fourier transform ( DFT).

CFT can be used in case of continuous aperiodic signals  while DFT for discrete aperiodic signals .

On the other hand, Laplace transform can be used in case of continuous signals and Z transform for discrete signals.

So I want to ask that can use of Fourier transform be minimized completely with the help of Laplace and Z transform?

The short answer is yes, if you have the Laplace or Z-transform of a function you do not need the Fourier transform.

This is because the CFT is a special case of the Laplace transform and the DTFT is a special case of the Z transform. The Fourier transform is used to find the complex sinusoids that compose a function, whereas the Laplace transform finds all complex exponentials that compose a function.

Knowing this, it is simple to derive the Fourier transform from the Laplace transform. The Laplace transform is defined as:

$$
\int_{-\infty}^{\infty}e^{-s}f(t)\,dt\\
s = \alpha+j\omega
$$

To find the Fourier transform from the Laplace transform, all you must do is assert $\alpha = 0$ for all $\omega$ and $t$. This evaluates the Laplace transform along the imaginary axis only, as any component lying directly on the imaginary axis is a complex sinusoid. This leaves us with:
$$
\int_{-\infty}^{\infty}e^{-j\omega}f(t)\,dt\\
$$

Which is in fact the Fourier transform.

A similar relationship applies between the DTFT and the Z-transform. The Z-transform is defined as:

$$
\sum_{n=-\infty}^\infty x[n]Z^{-n}
$$
Where $Z$ is a complex variable. To find all complex sinusoids in the discrete world we must constrain $Z$ to fall on the unit circle:
$$
Z = e^{j\omega}
$$

Which leaves us with

$$
\sum_{n=-\infty}^\infty x[n]e^{-j\omega n}
$$
Or the discrete time Fourier transform.

The reason the Fourier transform exists is because there are a lot of applications that only need the real frequencies of a function. If only the real frequencies are needed there is no point calculating the rest of the Laplace domain. Why do more work than necessary?

For more information check out the Wikipedia articles https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform#Relationship_to_the_Z-transform
and https://en.wikipedia.org/wiki/Laplace_transform#Fourier_transform And I would also strongly recommend the text book Linear Systems and Signals by B.P. Lathi for a much more indepth discussion on Laplace, Z, and Fourier transforms.


How to filter a signal knowing the frequency response of a system in MATLAB? I have a signal of period $T_0=8$, let's say $x(t)$, and it has the following Fourier coefficients:

$$ a_k=\frac{1}{4} \mathrm{sinc}^2\left(\frac{3k}{8}\right) e^{ik\frac{\pi}{2}} $$ for $1\le|k|\le6 $ and $a_k=0$ for the rest.

I also have a system which has a frequency response :
$$H(j\omega)=\frac{j\omega}{j\omega+\frac{3\pi}{2}}$$

Using MATLAB, what's the best way apply this filter to this signal?
I got the signal $x(t)$ from the coefficients and now I don't know how to proceed from here. I'll leave my code here:

T=8;
w0=2*pi/T;

t=0:0.001:2*T;
exp_jwt=exp(1j*w0*t);

k=-6:6;
ak=0.25*(sinc(3*k/8)).^2.*exp(1j*k*pi/2);
ak(k==0)=0;

x=zeros(size(t));
for (i=1:length(k))
    x=x+ak(i).*(exp_jwt).^(k(i));
end

Your signal is given by

$$x(t)=\sum_{1\le|k|\le 6}a_ke^{jk\omega_0t}\tag{1}$$

with $\omega_0=2\pi/T$.

An LTI system's response to $a_ke^{jk\omega_0t}$ simply is $H(jk\omega_0)a_ke^{jk\omega_0t}$, where $H(j\omega)$ is the system's frequency response. Consequently, the response to $x(t)$ is

$$y(t)=\sum_{1\le|k|\le 6}H(jk\omega_0)a_ke^{jk\omega_0t}\tag{2}$$

If you want to use Matlab, you could better use it to evaluate Eq. $(2)$, rather than to simulate the continuous-time system.


Difference b/w Spectral Entropy and Flatness Measure What is the difference between Spectral Entropy and Spectral Flatness?
As far as I understand, both the terms can be calculated from the power (or amp) spectrum, and describe the flatness (whiteness) of the spectrum.

If you wanted to distinguish the voiced activity from the unvoiced one using either of these measures, which one would you use and why?

Assuming that you calculate the Power Spectral Entropy as per the first link I've added to the question:

$${\rm PSE} = -\sum_{i=1}^N p_i\ln p_i $$

then the spectral flatness (using the same notation) would be

$$
{\rm Flatness} = \exp\left(\frac{1}{N}\sum_{i=1}^{N} \ln p_i\right)
$$
because the $p_i$ are normalized to unit sum, so the denominator in the second linked to question is one.

As you can see, the two are similar in form except that Flatness appears to a) exponentiated and b) assume that the $p_i$ are uniform.

As to which one to use, I'd probably opt for the PSE rather than Flatness. It appears to make more use of the spectral information.

To test my theory, I've written a comparison that starts off with a purely harmonic signal and goes to a purely white noise signal and calculates both measures for each blend.

There is very little difference between the two measures, as you can see from the closeness of the red and blue curves in the bottom plot.

Note that I've normalize both curves so they go between 0 and 1, just to make it easier to compare them.

Bottom line: Choose the simplest; there appears to be little difference between the two for this example.

R Code Below

#30534
T <- 1000
white_noise <- rnorm(T,0,1)

K <- array(seq(1,5),c(5,1)) %*% array(1,c(1,T))
omega <- 2*pi*0.053492384
t <- array(seq(0,T-1), dim=c(1,T))
t <- array(1,c(5,1)) %*% t
harmonic <- colSums(sin(omega*K*t))

spectral_entropy <- function(signal)
{
  psd <- abs(fft(signal))^2
  psd <- psd/sum(psd)

  pse <- 0
  for (ix in seq(1,length(psd)))
  {
    pse <- pse - psd[ix]*log(psd[ix])
  }

  return(pse)
}

flatness <- function(signal)
{
  psd <- abs(fft(signal))^2
  psd <- psd/sum(psd)

  flatness <- 0
  for (ix in seq(1,length(psd)))
  {
    flatness <- flatness + log(psd[ix])/length(psd)
  }
  return (exp(flatness))
}

P <- 100
pse <- rep(0,P)
flt <- rep(0,P)
for (alpha in seq(0,P))
{
  signal <- white_noise*(alpha/P) + (1-alpha/P)*harmonic
  pse[alpha+1] <- spectral_entropy(signal)
  flt[alpha+1] <- flatness(signal)
}

normalize <- function(signal)
{
 return( (signal - min(signal))/(max(signal) - min(signal)))
}

layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
plot(abs(fft(harmonic)), type="l", col="blue")
title("Harmonic Spectrum")
plot(abs(fft(white_noise)), type="l", col="blue")
title("White Noise Spectrum")

plot(seq(0,P)/P, normalize(pse), type="l", col="blue")
lines(seq(0,P)/P, normalize(flt), type="l", col="red")
title('Comparison of PSE vs Flatness')


How can I duplicate this type of image? I'm very new to the idea of signal processing, could someone point me in the direction of of a program to toy with images, below is a link to 2 taken a split second apart that showed up on my Amazon app. I've been making decimated images recently in Photoshop and want to play with images in a way that may have more effect like this

The second image is sort of pixelated but not totally decimated image of the first image, an almost plain book cover. It's a split second of the lower one of this pixelation  then snaps into the regular image. How can I take original images and process them with different signals to this effect. Notice there's banding on the top and bottom so it's not just a simple rescaling in Photoshop which I'm very familiar with.

I just cropped the original gray image:

and converted it to JPEG using irfanview, settings to 15% (my first guess, good feeling isn't?):

You can see see the same yellowish shades on the right, due to heavy chroma quantization. Specifying that the image should remain grayscale gives you this:


Objective, repeatable sound quality measurements and metrics I am currently doing an electronics design for a basic 2.5W speaker amplifier for an embedded system based on COTS parts (e.g. codec and amp), similar to what one would find in a tablet.

Aside from a wide range of parts to choose from, there are other engineering trade-offs, such as the filter on the CLass D amp for EMI/EMC and the noise on the power supply rails.

This, of course, is even before deciding on a speaker.

Is there a standard method for an objective and repeatable evaluation process? If so, what sort of equipment and methods do I need, and what metrics do I need to calculate, and how to these metrics relate to the perception of audio quality?

(Yes, of course, there is SNR THD and frequency response, but I am yet to relate these to how a human perceives the quality of the audio).

(Second note: This isn't an evaluation of codecs, compression, sampling or other DSP artefacts.  I have control of the bitstream I send to the audio codec chip. Advice of choosing the signal I should send is most welcome).

i dunno anything about EAQUAL or PEAQ.

judging which sounds better is quite subjective, so much so that i consider it nearly a useless question.

but if you have a reference (that you define as "good"), and you have several options to compare to the reference, and your goal is find which option most sounds like the reference, then i would recommend AB testing (a little bit different from ABX testing).  at least blind, and better yet, double blind.

the subject listens to A and to B.  as many times as they would like.  pressing A always gets them sound A and pressing B always gets them sound B.  then the only question to ask the subject is "Do these sound the same or are they different?"

put in an equal number of pairs that are truly identical as those that are the difference you're trying to test.  so for two different sounds, the subject hears AA, AB, BA, and BB.  if the subject scores more true positives than false positives and more true negatives than false negatives, he or she might very well be hearing the difference.  get lotsa subjects and have them hear lotsa different sound pairs and you can get some meaningful results.


Do I need the sampling frequency in order to gain periods out of FFT? Currently I am using the transformInPlace() method from here.

for a project at university.
First of all I thought FFT would be the perfect solution for my problem, but
now I am not so sure anymore.
I have to alter a certain component of a software called LIMBO.
LIMBO is a load intensity modeling tool.
Every webpage has a trace that records the number of site views during a fixed period of time.
This leads to trace files like this data.

You can see lots of (time, arrivalRate(time)) pairs. This is all LIMBO gets as input.
The goal is to search for dominant periods within the noisy trace via the FFT.
The problem are the time values of the (time, arrivalRate(time)) pairs.
For a fixed time (for example 1 day) the page counts views. The total number of views than is
stored inside arrivalRate(time).
Then the page listens for another day in order to calculate the next arrivalRate(). This process is repeated
the whole time.
So, the time values inside the trace (for example 0.5,1.5,2.5) are not necessarily days (They  could be seconds, weeks or even months as well). They are abstract and depend on individual fixed listen period of every webpage mentioned above. (probably listen period  is called sampling frequency but I am not sure)
Therefore the information about the listen period of the page can not be retrieved out of the trace file.

My question is:
  Is it a problem to know nothing about the listen period of the given page if I want to calculate and interpret dominant periods
after a FFT has been applied on the trace file?

For example I use this formula to calculate $f_i $(the frequency of the $i$-th array position)

$i$: current position in the array

$N$: length of the array

$f_s$: the sampling frequency that is probably not revealed in the trace (so far I just used values like "1" or the median of all arrival rates in the trace)

$i = N f_i/f_s$   (so I use $1/f_i = {\rm period}_i = N/(f_s i)$ for calculating the $i$-th period)

This is DSP – everything that defines time or frequency is only measurable with respect to your sample period; so without that, your FFT output is but a relative measurement.

In short: If you can live with giving information like "events per sample period", then you're fine, if you need to say something about "per hour" or so, well, that information has to come from somewhere, and if you're not feeding in that information into your system, it can't have it.


Detecting meandering rivers in images I have $n$ surfaces: $z_i(x,y)$ with a measured attribute (variable) on each surface: $a_i(x,y)$. Most of the surfaces will have a random distribution of the attribute across the surface, but some surfaces (the interesting ones) will show a meandering river pattern:

I need your help in coming up with a measure that will tell us which of the $n$ surfaces are most likely to have such a pattern.

There are many possible maps with the same histogram as shown below; so the measure needs to "reward" spatial continuity.
To illustrate this I have created a random image with nearly the same histogram as the river image:

So image statistics ala entropy may only be part of the solution.

Here is an example of an image without a meandering river pattern:

My images are synthetic (made in Matlab). In real life the image without the pattern may have somewhat more spatial continuity in the form of small blobs of similar value.

Here are the images in grayscale:

A very simple measure would be to compare each row in the image with the row above it, allowing some horizontal shift.

I've hacked together this simple algorithm in Mathematica:

Mean[MapThread[
  Function[{line1, line2},
   Min[Table[Norm[line1 - RotateLeft[line2, shift]], {shift, -5, 5}]]
   ], {s[[2 ;;]], s[[;; -2]]}]]

It simply takes each pair of adjacent rows, rotates one of the rows by -5..5 pixels, and takes the smallest euclidean distance. This yields one euclidean distance for each row pair. I simply take the mean (but depending on your actual data, a truncated mean or median might be more robust).

These are the results I get for artificially generated samples (Formula: Normalize(random noise * (1-factor) + signal * factor))

If I plot the result against the signal strength, the algorithm seems to measure the "meandering river signal strength" quite well:

EDIT: I forgot to normalize the input samples. Fixed that an uploaded new result images


FFT calculation of an N-tap FIR filter I have the following question:

a) For an N-tap FIR filter, sketch the direct form realization. How many multipliers per sample? You may choose 8.
b) For an FFT calculation of an N-tap filter sketch the realization. (N=8, FFT size is 256)

The answer to a) is given here, and I think 8 multipliers are needed per sample, assuming an 8-tap filter.
With respect to b) I am confused. I don't know how to form a realization of the FFT of a filter. Could somebody help me out with b). Thanks!

hint only, since this is homework and you could get a lot out of this!
The mathematical operation a filter performs between filter impulse response and signal is called __________(1). When we process both filter and signal with an FFT (which is just a fast implementation of __________(2)), then that operation on the originals is equivalent to __________(3) after the (2).
To do (3) on both (2)'ed vectors, these vectors need to be of the same  __________(4), here 256. Often, the signal itself is longer than 256 - (length of filter). This is a problem, because (1) under the (2) is not  __________(5) as the filtering operation normally is, but  __________(6) due to the periodic nature of the (2).
To alleviate this, a segmented approach is chosen.


Is there a difference between Image Alignment and Stereo Rectification for stereo correspondence or matching? As far as I know,
Image alignment just finds the best fit given two images of the same scene.

Whereas Image rectification, with regards to stereo vision, warps the images so that the epipolar lines are all on the same level. Essentially all common pixels between the two images are warped to be on the same y axis.

My question is, what is the difference between Image Alignment and Image Rectification, given a pair of images, when doing stereo correspondence?

If I were to find the disparity between two images, how is image alignment different to image rectification?

There two ways to look at this problem.

In simple terms, image rectification warps both images onto a common coordinate frame by typically estimating the transformation using the epipolar geometry.
Image alignment finds the transformation from one image to the other. It doesn't guarantee any constrains on the epipolar geometry and only one single image is warped.

Image rectification might benefit from image alignment (of local patches) to obtain the correspondences. Once these are found, the epipolar lines can be forced to be parallel.


channel estimation based on time-domain pillots I am confused if it's possible to estimate the channel based on the time-domain pilots. I mean I have such signal $x$ with length $N$ after taking its $iFFT$ we get $X = ifft(x)$. I have perfect knowledge on pilots signal in time domain, I mean I know what is $X_{1:4:N}$ (Again, it's in time domain after taking the ifft). The signal $y = h*X$ is the received signal representing the convolution between signal $X$ and channel $h$. Taking the $FFt$ to signal $y$ resulting $Y = fft(y) = fft(X)*fft(h)$.
My issue, if it's possible to estimate the channel in this case based on the pilots known in time-domain which are $X_{1:4:N}$ ?

if we ignored the OFDM system for the moment, can we estimate the channel h based on y in my above mentioned system?

Not really, because what you're describing is like if you could sample the channel, but at (symbol rate / 4), if I understood your approach correctly.
Hence, you get aliasing in your channel estimate, if it changes faster than (symbol rate / 4).
Your system, however, sees a channel that has full symbol rate (at least!) in bandwidth.
If your channel is slow enough, i.e. it's coherence time $T_C > 4T_S$, then this subsampling can be OK.
But you haven't given us any info that says it is (if you really ask us to ignore the OFDMness of it all).
If, on the other hand, we don't ignore the fact that this is supposed to be an OFDM system (in which case you really mustn't insert pilots in time domain, that breaks OFDM completely!), then your pilots are far too closely spaced, assuming your OFDM system has more than $N=4$ subcarriers: To use OFDM, your subchannels need to be coherent for one OFDM symbol duration $T_{Sym}=NT_\text{sample}=\frac{N}{f_\text{sample}}$.
I think you might be confusing something. OFDM systems of course use pilots on their subcarriers, which is in frequency domain, that they replace data symbols with every couple of OFDM symbols. But that means the pilots need to be inserted into your data stream, not into your time domain signal.


Difference between these frequency response curves I am thinking about a simple mass-damper system response. Most single DOF systems show a frequency response plot that looks like this:

Note how increasing the damping ratio always reduces the output response, and at high frequencies, the response approaches the same line regardless of damping ratio.

But, the way I am used to thinking about it, the response should really look like this:

Note how the increasing the damping ratio reduces the response below $\omega/\omega_0 \cdot  \sqrt{2}$, but increases at frequencies above $\omega/\omega_0 \cdot  \sqrt{2}$.
The response does not approach a downward sloping line as damping ratio increases.  Rather, it approaches $1$ as damping ratio increases.

So what am I missing here?  What is the difference in these two curves?

Thanks

The two frequency response plots are the result of two different models. In the first one, the model is a mass-spring-damper system as shown here, and the plot shows the ratio of the amplitudes of mass displacement and driving force.

The second plot shows the transmissibility of a single DOF mechanical system as shown below:

(from Fundamentals of Mechanical Vibrations by S.G.Kelly).

Note that here, unlike in the other model, the support is not fixed. The transmissibility shown in your second plot is the ratio of the amplitudes of the mass displacement $x(t)$ and the (forced) displacement of the support $y(t)$.

In formulas, your first plot shows the function

$$M(\omega)=\frac{1}{\sqrt{\left(1-\left(\frac{\omega}{\omega_0}\right)^2\right)^2+\left(2\zeta\frac{\omega}{\omega_0}\right)^2}}\tag{1}$$

where $\omega_0$ is the system's natural frequency, and $\zeta$ is the damping ratio. In electrical engineering terms, this is the magnitude of the frequency response of a second-order low pass filter.

The second plot shows the transmissibility given by Eq. (3.73) in Kelly's book:

$$T(\omega)=\sqrt{\frac{1+\left(2\zeta\frac{\omega}{\omega_0}\right)^2}{{\left(1-\left(\frac{\omega}{\omega_0}\right)^2\right)^2+\left(2\zeta\frac{\omega}{\omega_0}\right)^2}}}\tag{2}$$

Note that $(1)$ and $(2)$ have the same denominator but different numerators.


Simulate time series given temporal auto-correlation functions Given a random process $x[n] \in \mathbb{R}$ (say of length $N$) and all correlation functions such as:

\begin{align}
\langle x[i]\rangle\\
\langle x[i]x[j]\rangle\\
\langle x[i]x[j]x[k]\rangle\\
\vdots
\end{align}

Is it possible to simulate a single realization of the random process and if so how?
Thinking abstractly about the problem it seems like any trajectory could be mapped onto a vector in $\mathbb{R}^n$ and there is some probability distribution on $\mathbb{R}^n$. Sampling this distribution and then picking out the corresponding trajectory would then constitute of "simulation" of a single realization that I am looking for. It seems like given information about all of the auto-correlation functions it should be possible to determine the probability distribution. This seems related to the inverse moment problem.
What if we truncate and only have up to $k$-time correlation function with $k<N$? My guess is that it is not possible to uniquely generate the probability distribution but it should be possible to (non-uniquely) simulate a process with matching auto-correlation functions up to $k$.

Here is one possible method that you could possibly use.

Let's model your random process $x[n]$ as an auto-regressive (AR) process of order $P$, that is,
$$x[n] + \sum_{p=1}^P a_px[n-p] = c + \epsilon[n]$$
where $c$ is some constant that I will set to 0 in what follows for the sake of simplicity, $\{a_p\}_{p=1}^P$ are the parameters of the AR model and $\epsilon[n]$ is some white noise of variance $\sigma^2_\epsilon$.

Multiplying each side of the above equation by $x[n-l]$ for some $l \ge 0$ leads us to
$$\left(x[n] + \sum_{p=1}^P a_px[n-p]\right)x[n-l] = \epsilon[n]x[n-l].$$
Now, taking the expectation on each side, we obtain the Yule-Walker equations
$$C_x[l] + \sum_{p=1}^P a_pC_x[l-p] =
\begin{cases}
0 & \text{if} & l > 0, \\
\sigma^2_\epsilon & \text{if} & l = 0
\end{cases}$$
using the fact that $\epsilon[n]$ and $x[n-l]$ are uncorrelated for $l > 0$ (you can see that from the equation of the AR process above). Using this equation with $l = 0, 1, 2, \dots, P$, you have $P+1$ equations to solve for the $P + 1$ unknowns (i.e., $\{a_p\}_{p=1}^P$ and $\sigma^2_\epsilon$).

You then have a model to generate a time series given the temporal auto-correlation.


You know how to map from dBm to mW and the other way around so it seems to me this is a lossless transformation and you should be good to go by converting the dBm value to the mW value and doing the distribution fitting. You can then use it for whatever you want but just remember that you have done the conversion so will need to convert back to dB scale if that is what you want to report in the end.

Edit:
From the comment, since MATLAB is complaining about too small of values I suggest just converting to microwatts. So do:

z = 1000*z;
pd = fitdist(z, 'Rician');

If it still complained, we could just keep moving down to nanowatt and so on. As long as you know that when you use that distribution you need to keep the units straight, then there is no harm in going from dBm to mW, $\mu$W, nW, and so on.


MATLAB : How to filter a data signal I have a data signal (column vector) with N values taken over at specific sample rate (dx) at a fixed frequency. After taking the Fourrier transform of this signal and representing it in the Fourier Domain (with fftshift), I need to eliminate the parts situated in the origin and in the positive side. The final signal I need is the one with the negative shift. Therefore, after filtering the unwanted parts I will need to perform an inverse Fourier transform to recover the signal in time domain and center it in the origin.
Finally, I do not understand how to filter the spectrum when there is no frequency variation and when I’ve tried specifying a sample frequency and then applying a matlab filter I didn’t get any results. Can anyone help me to understand how to make this work please?
PS: the spectrum is in the figure below

This task its not so hard, you only should create a vector, using ones, of the same size and set to zero the positions that correspond with the frequencies that you don't want. Apply an element-element multiplication between the vector and your signal. Finally apply the inverse Fourier transform in the obtained vector. The response you will obtain will be a complex one, so you need to graphic only real or imaginary part if you want to appreciate it in a graphical way, you also can graphic the absolute of the signal in order to observe the magnitude of it.

This type of filter is known as ideal filter and is impossible to apply in real world. In Fourier domain, filter is just a multiplication while in space domain to filter is necessary to apply a convolution between the signal and a kernel.


Change pitch of mp3 audio on a microcontroller I have an MP3 decoder running on a ESP32, and I'd like to change the pitch of the streaming audio in real time.

I plan to build a little toy that streams Amazon Alexa and changes her voice so that it is higher pitched. So no pro audio requirements at all, it is just supposed to be a funny effect.

I was thinking maybe I can operate on the raw DCT data? Or is there another way that is fast enough for that chip? I cannot use time compression/dilation because the audio is streamed.

I'm using the MAD fixed point decoder.

A description of (or a link to) an algorithm working on DCT or PCM data would be fine, but code or pseudo-code in whatever language would be very welcome.

Use a "phase vocoder" which will use short term FFT properties to change the pitch without changing the duration.

There is plenty of information and links to "phase vocoder" on Wikipedia and Google in general.

It is a great educational exercise to go through its operation as you learn from it the critical importance of the phase information in the FFT for proper reconstruction.


Energy calculation in frequency domain I was just wondering... The formula I learned to calculate the energy of the signal is expressed in the time domain:

$$E_x^{\text{time}} = \sum_{n=-\infty}^{\infty} |x[n]|^2$$

Then, what does the amount of energy gotten from the magnitude spectrum mean?

$$E_x^{\text{frequency}} = \sum_{f=0}^{fs/2} |X[f]|^2$$

Supposing that both amplitudes are expressed in decibels, $fs/2$ is the maximum representable frequency (by Nyquist theorem), and $E_x^{\text{frequency}}$ is taken over the same time frame as $E_x^{\text{time}}$. Should the quantities be the same?

You need Parseval's theorem. For the discrete-time Fourier transform (DTFT) you have the following relation:

$$\sum_{n=-\infty}^{\infty}|x[n]|^2=\int_{-1/2}^{1/2}|X(f)|^2df\tag{1}$$

where $f$ is normalized by the sampling frequency. For the DFT you have

$$\sum_{n=0}^{N-1}|x[n]|^2=\frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2$$

So due to Parseval's theorem it is always possible to compute a signal's energy in the time domain as well as in the frequency domain.


Cross Correlation dimensions I have been trying to do the cross correlation of two signals, and I am expecting the output to have the same dimensions as the two signals have. However what I am seeing is that the cross correlation has twice the amount of points

The red here is the cross correlation of the green and blue. The two starting signals cover $1$ hour of time, sampled at $100Hz$, where as the cross correlation seems to represent $2$ hours, what I want to see is how correlated they are over the 1 hour. Can anyone explain this? (also apologies the graphs look like the joker)

When you compute cross correlations, you shift one signal over the other and compute the (normalized) inner product of the overlapping sections of the two signals. The point is that you can shift in both directions, i.e., your time shift (in samples) can be positive or negative.

Assume both signals have $N$ samples. Your central data point in the cross correlation is when you haven't shifted your second signal at all, i.e., all $N$ points are overlapping. Now, you can shift "to the right" by 1 sample: your overlap is $N-1$. Shift another one: your overlap is $N-2$. Keep going, until you shifted $N-1$ samples: there is only an overlap of 1 sample. This gives one half of your correlation function. You get the other half by shifting in the other direction, "to the left", you get another $N-1$ samples like that. In total, this gives $2N-1$ samples, which explains the "doubling" you observe.

Example: Correlate the sequences [1,2,3] and [4,-5,6]. Normalization left aside, we get as inner products:

[1,2,3]'*[0,0,4] = 12
[1,2,3]'*[0,4,-5] = -7
[1,2,3]'*[4,-5,6] = 12
[1,2,3]'*[-5,6,0] = 7
[1,2,3]'*[6,0,0] = 6

Therefore, from N=3 original lags, we get 2N-1=5 lags in the correlation function.

If you have a good reason to believe that your signals are periodic, you can use a cyclic convolution, which will result in N samples exactly. But his only makes sense for periodic signals.


Why look at power spectral density for stochastic processes? I have been told that for deterministic signals, it makes sense to look at their respective Fourier transforms/spectra.

For stochastic processes on the other hand, I am supposed to work with power spectral density in terms of qualitative analysis.

Why?

Because a stochastic process itself doesn't have a Fourier transform.

That's really all there is to it.

You can only transform signals (i.e. functions over a body isomorphic to $\mathbb R$, for example, functions of time). You can't transform a random variable whose individual realizations are such functions!


How would I represent an image using a basis $A$, given a sparse or compressible image? Looking for a Basis representation:
Given a sparse or compressible image M, how would one represent this image using a basis A?

If I let m = vec(M)

m = Ax

where A is my "representation" basis and x is the vector of coefficients in the A domain, what is the method for actually getting A?

The textbook algorithm is to minimize

$$\|Ax-m\|_1+\alpha\|x\|_0$$

where $\|x\|_0$ counts the number of non-zero coefficients. Since that is a bad function for the determination of a descent direction, one softens the counting norm to

$$\|x\|_ε=\sum_k|x_k|^ε$$

for some small value of $ε$. For instance, for $|x_k|>\exp(-\tfrac1{\sqrtε})$, which is a very small threshold, one gets $|x_k|^ε>\exp(-\sqrtε)\ge 1-\sqrtε$, so that all not too small and not too large values get mapped to $1$ and only really small values get mapped to zero.

To get a really smooth function, mostly for theoretical purposes, one can smoothen the kink at $0$ by using

$$\|x\|'_ε=\sum_k(\delta^2+|x_k|^2)^{ε/2}$$

and obviously choosing $δ$ small enough so that $x_k=0$ still gets mapped close to $0$, for instance by choosing $δ=\exp(-\tfrac1{\sqrtε^3})$.

Now that this answer is completely orthogonal to the question in its current form, I hope that someone else has more insight into the topic of matrix selection. What I know is this:

In the beginning, before compressed or compressive sensing existed as word, $A$ was a matrix of wavelet or wavelet-like coefficients. The whole idea of wavelet based compression (and JPEG) is to zero out small coefficients (and truncate the binary representation of non-null coefficients). This can be improved upon in one direction by observing and utilizing correlations in the detail coefficients over multiple scales, or in the other direction by making the basis vectors even more uncorrelated and thus reducing those fractal features.

There may also be algorithms that start with a random matrix and modify or reject it based on certain frame conditions.


number of subcarriers for an OFDM system For an OFDM system, it is necessary to know the number of subcarriers. For example for the subcarrier $k^e$, we must have: $B_k <B_c$ with $B_c$ the channel coherence bandwith. And $B_k=\frac{1} {NT+I}$ with $I$ is the guard interval, $T$ is the symbol time and $N$ the number of subcarriers.

Here is my questions:

For example, for $N=3$, we have $B_k<B_c$. Is the number of subcarriers small?  Because I always see, N = 64, 1024, ...
Is there a minimum number of subcarriers?

You've got that right – when designing an OFDM system, you pick the number of subcarriers usually just large enough to make the subchannels look flat.

You avoid overly large subcarrier numbers, because more subcarriers means longer symbols, and you need to restrict the symbol length to stay below coherence time of the channels.

Now, 3 is indeed a choice that is very small. I wouldn't build an OFDM system with 3 subcarrieres, simply because the smallest cyclic prefix that you can build is then $\frac 13$ of the symbol length – and chances are, you need more than one $\frac1{B_{subcarrier}}$ as guard interval to "capture" the whole channel impulse response.

Then again, your channel seems to be very benevolent, since it has such a high coherence bandwidth, so maybe that's really the case.

However, in a situation where you'd do 3- or 4-subcarrier OFDM, you'd simply not do OFDM alltogether and would probably either simply use a single-carrier system with a proper equalizer (might be easier, even), which will not waste as much spectrum as OFDM, or you'd go for filterbank multicarrier systems, because you will have a very bad time¹ trying to achieve synchronicity in such an OFDM system. Schmidl&Cox will be very ambiguous prior to first lock, and the SNR of any timing symbols will be less-than-good with just three values to work with per OFDM symbol.

¹ yes, that's a pun.


IIR Filter for Smoothing (Low Pass Filter) I am using IIR filter for smoothing

$$y[n] = ax[n]+(1-a)y[n-1]$$

My question is, if I add another IIR filter, will it be the second order of IIR filter? If not, what it can be called?

My second filter is

$$y_2[n] = ay[n] + (1-a)y_2[n-1] $$

If you apply two filters in a series cascade, then the behavior of the cascade can be expressed in two different ways. In the time domain, the overall system's impulse response can be calculated by convolving the impulse responses of $y[n]$ and $y_2[n]$ together. For IIR filters, this can be somewhat cumbersome.

In the frequency domain, the overall system's $z$-domain transfer function can be calculated by multiplying the transfer functions $H_y(z)$ and $H_{y_2}(z)$ together. This is usually a much easier route for filters with feedback.

In your case, the two filters actually have the same input/output relationship (assuming that $y[n]$ is the input to $y_2[n]$. Using the $z$-transform, it's easy to find that:

$$
H_y(z) = H_{y_2}(z) = \frac{Y(z)}{X(z)} = \frac{a}{1 - (1-a)z^{-1}}
$$

Using the relationship I mentioned above, you can calculate the transfer function of the two filters in cascade using:

$$
H(z) = H_y(z) H_{y_2}(z) = \left(\frac{a}{1 - (1-a)z^{-1}}\right)^2
$$

$$
H(z) = \frac{a^2}{1 - 2(1-a)z^{-1} + (1-a)^2z^{-2}}
$$

We can just as easily use the inverse $z$-transform to arrive back at the difference equation for the two cascaded filters:

$$
y_c[n] = a^2x[n] - 2(1-a)y[n-1]+(1-a)^2y[n-2]
$$

By inspection, we can state that this is a second-order filter (provided $a \not= 1\ $), as you suspected.


Generating a Sine Signal (for example a voltage sine wave) with changing frequency and changing sample rate (time)? i have a measured signal (e.g. voltage sine wave) and i would like to generate a sine wave which matches this meausred signal (with the generated sine wave i remove the harmonics from my measured signal), just like a PLL block in simulink.

for that i have to do the 2 steps:

1- Calculate the frequency of the measured signal from detecting the Zerocrossing.

2- calcualte the sine function using sin(2*pi*f(t)*t).

I have the following problem.

1- calcuating the frequency by counting the samples between two zerocrossing would results in getting the frequency at one point so my calcualted fequency does not match my sampled measured signal.

2- even when interpolating or using the repmat function i am able to get the f(t) to be the same length as the time (t) vector yet my sine wave is not correctly calcualted. i know that:

I am adding the zerocrossing times on my time vector which makes the last one inconsistance which mean the sampling rate is changing also.

so how to genereate a sine wave which changing frequency and changing sampling rate of the time signal?

Best Reagrds

If you are doing this in realtime, then your generated signal WILL ALWAYS lag behind the original.
Your sampling rate will not change because of anything your calculation does.
You have a problem with discontinuities when the frequency of your measured signal changes - you are abruptly changing from the old to the new frequency.  This will cause skips in the generated signal.
You will not always be able to exactly match the real signal - the zero crossing may take place between two samples.  This will cause your calculated frequency to be wrong.
The calculated frequency won't necessarily fit well with your sampling rate, so that it won't end exactly on a sample - this will cause another skip in your generated signal.

You will have to gradually change f(t) in order to avoid the abrupt changes.  So, you measure f one time and it is (say) 1000Hz, then the next time it is 1100hz.  You cannot simply change from one to the other.  You will have to "walk" the value for f from 1000Hz to 1100Hz over the length of time it takes for one half cycle of 1100Hz.

That should improve things a good bit.


Need a better step detection algorithm I have a time series with lots of steps/jumps (data file here). A plot is given below. I would like to subtract an appropriate value for each of these square wave features to bring them back down to the baseline of the signal. A median filter works really well for removing a small number of outliers in a row, but in this case I probably need a different approach since the square wave jumps can have different durations as seen. A common method I've seen for doing this is to compute first differences of adjacent samples, and look for large differences to detect jumps. I implemented this method but the problem is it often fails, since the one tunable parameter for the method is a threshold value $t$ which the first differences must cross in order to detect a jump:
$$
| x_{i+1} - x_i | > t
$$
As can be seen in the plot below, the jumps I have are often different sizes, so a constant threshold value isn't the best approach. In particular, in some cases there is an interesting signal where adjacent samples can change by large values without being a jump! I have highlighted such a region in red.

Below is a zoomed in view of the red box area. You can see there is a square wave jump followed by an interesting signal. The red arrows depict a place where adjacent samples from an interesting signal have a larger distance between them than some of the jumps in the signal. Therefore a constant threshold method with finite differencing will not work for me.

Does anyone know of a robust procedure to detect and subtract the square wave jumps to end up with a smoothly varying signal with no jumps? I'm sure this must be a solved problem but I haven't had much luck searching online.


The Fourier Series Of This Triangle Wave I am using matlab to study digital signalling and have come across a problem which i was wondering if anyone with more experience could help me with.

I need to work derive the Fourier series of a triangle wave that i have generated, I just do not know how to actually go about this problem in Matlab.

I am generating a 100hz Triangle signal using the following code:

t = 0:1/10000:1;
f=100;
x1 = sawtooth(2*pi*f*t, 0.5);
plot(t,x1);
axis([0 0.10 -1 1]);

Now how should i go about deriving the Fourier series of this signal, i am completely lost.

Any help would be appreciated.

Check out Wikipedia's Fourier Series page.  Their "Example 1" shows how to derive the Fourier series of a sawtooth wave.  All you have to do is normalize the results for your particular time and amplitude values.


Image Segmentation - Comparing GrabCut Method vs. Level Set Method Checking the literature about image segmentation, if I am not wrong, there are two main approaches:

Variational approach e.g. level-set optimization in Chan-Vese
Combinatorial approach e.g. using graph cut in GrabCut

Can someone give me some high-level comparison between these two directions?

For example, I feel with level-set optimization one can easily  incorporate image structure (e.g. edges) into the energy function, which is not easy in Grabcut. Am I correct?

Thanks.

Level set is for continuous segmentation, multiphase. But Graph cut is for discrete segmentation, only two phase


Why do I have to divide by $\pi$ in order to scale amplitude of sine frequency spectrum? Recently I found my old code to transform sine wave from time domain to frequency domain. Below I post simplified version of it.

from math import sin, pi
from matplotlib.pyplot import *
from pylab import zeros, fft

fs = 500
duration = 2 * pi
amplitude = 2
frequency = 25
values_in_time_domain = []

for t in range(int(duration * fs)):
    values_in_time_domain.append(amplitude * sin(frequency * t / fs))

subplot(2, 1, 1)
title('Sinus from 0 to 2 pi in time domain; A = 2, f = 25, fs = 500')
plot(values_in_time_domain)

transform = fft(values_in_time_domain)
transform_abs = [abs(x) for x in transform]
values_in_freq_domain = [transform_abs[f] / fs / pi for f in range(int(fs / 2)+1)]

subplot(2, 1, 2)
title('Above sinus in frequency domain')
plot(values_in_freq_domain)

show()

Why do I have to divide every element from transform_abs by pi in the following: transform_abs[f] / fs / pi?
I noticed I need to do it in order to have properly scaled $O_y$ axis.

The FFT should be divided by the number of elements in the vector.  You specified the vector length to be duration = int(2*pi*fs) which is a really weird way to defining the length, btw.  Thus, you need to divide by int(2*pi*fs).  Yes, you need to divide by a 2 as well, because the signal energy is distributed evenly between positive and negative frequencies.


PSS convertion from Frequency to time domain Hi I have a question regarding primary sync signal in LTE.

It is generated in frequency domain and its length is 63. How to convert this signals frim frequency to time domain. Should I use IFFT? what is the length of IFFT and the signal in time domain?

-Thanks

The PSS signal is generated based on Zadoff-Chu sequence and is mapped to the central 62 sub-carriers apart from DC sub-carrier as explained here. If you want to view the PSS alone, leaving aside other subcarriers in the PSS containing OFDM symbols then yes, you'll have to do a 64 IFFT on the length 63 Zadoff-Chu + DC symbols. It is worth nothing that on the receiver side, synchronization is achieved by locking on to the central 63 subcarriers using a 64 point FFT.


Power/Amplitude Response of a system I am currently looking over filters and understand the frequency response of a system, e.g. a low-pass filter. However I am confused by the 'Power/Amplitude response' of a system as detailed in my notes below:

  Square of the amplitude response, or power response:

  $$|H(\omega)|^2=H(\omega)H^*(\omega)=H(\omega)H(-\omega) $$

I created a basic matlab scenario for the a low-pass filter and plotted the frequency response and the power response:

From the top graph I understand that as the frequency increases along the x-axis you get an attenuation for difference frequency components in db.
However I do not understand what the significance or intuitive value is gained by squaring the attenuation (or gain). I can see that the gain becomes more negative as frequency increases and the amplitude response increases.
So for example my 3dB cut-off point is at 1rad/s, 10^(0). This has a corresponding amplitude response value of 9. Does this mean the amplitude has decreased by a ninth? Because surely if my gain were positive my amplitude response value would also be a positive number and so the opposite cannot be true?

Any help would be greatly appreciated!

Regards

Pat

Your second plot is the square of the magnitude of the frequency response in dB, which indeed doesn't make much sense. You must use $|H(j\omega)|^2$, and not $(20\log_{10}|H(j\omega)|)^2$, as you did. Of course, $|H(j\omega)|^2$ will not give you any new information compared to $H(j\omega)$.


In block codes (channel coding), in practice how is the message length ($k$) fixed? In block codes (channel coding), in practice, how is the message length ($k$) fixed?  If $k$ is large, the practical implementation of it requires huge memory if look up table approach is adopted.

What is the disadvantage with short $k$'s?  Code rate will degrade?

There are a few rules of thumb:

For a fixed code rate $r = \frac{k}{n}$, the amount of gain you observe from applying coding typically increases as you increase the block size $n$. Accordingly, the complexity of decoding the larger blocks is often larger.
If you decrease the number of information symbols in the block $k$ and keep $n$ fixed, the code rate decreases. This will result in higher coding gain, but lower information throughput for a given channel symbol rate.
If you increase the number of information symbols in the block $k$ and keep $n$ fixed, the code rate increases. This will result in lower coding gain, but higher information throughput for a given channel symbol rate.


Nyquist Frequency Confusion 1- If I have a sine wave with period of $T$, I need to sample at least every $T/2$ to be able to reconstruct the sine wave. Let's look at this:

This way I'd get a series of $0$s and all information seems to be lost. Where am I wrong?
$x(t)=\sin(2\pi ft)$

Let $f_s=2f$; then

$$x[n]=\sin(2n\pi f/2f) = \sin(n\pi)=0.$$

2- Assume I have sine wave of frequency say $100\text{ Hz}$ , and I sample it at $101\text{ Hz}$, this would cause aliasing as we know. My question is, what happens in the frequency domain?

For the number of samples given, your initial statement is false.  Sampling at exactly T/2 is a mathematical limit.

For any finite number of samples, you need to sample a band-limited pure sinusoid of period T more frequently than T/2 to avoid aliasing or other pathologies, or at greater than 2X the highest spectrum frequency (for a baseband range). (and in real-life, to take into account the transition bandwidth and stop-band of any low-pass anti-aliasing filtering).  The shorter the sampling interval, the higher the sample rate needs to be above 2*Fmax (for a given baseband band-limited range) to reduce numerical, quantization, and other signal-to-noise induced errors to below some given threshold.

As you approach an infinite number of samples (beyond the life of the universe), you approach being able to reconstruct a sinewave with sampling at T/2.


In computed tomography (CT), why is 'Inverse Problem of Radon transform' studied? As I know, the Radon transform is a very important tool in CT by Beer's law. Thus, finding $f(x)$ of Radon transform $Rf(L):=\int_{L} f(x)dl(x)$ is helpful in CT. Nowaday, the Filtered back-projection formula (FBP) for all dimensions $n$ (which can find the $f(x)$) is obtained and is used already.
My question is here. Why is people studying 'Inverse Problem' in Radon transform although FBP is obtained?

Note that the original mathematical back-projection method assumes that the value of the Radon transform is known for all lines. This is an infinite amount of information, so it is pure mathematics at that level.
Since practical tomography works with a finite number of values of the Radon transform, the inversion is ill-posed. This cannot be ignored due to the finite precision of any numerical measurement and to the impact of noise.
In addition, since exposure to X-rays is harmful, it is imperative to continue to seek and develop algorithms that achieve results with fewer X-ray exposures.


How to Find the Frequency Response of a Communication Channel from Input and Output Symbols in MATLAB? I want to find the channel frequency response of a digital communications system. I have the functions of the input symbol (a triangle) and the output symbol - a distorted triangle: $$\dfrac{1}{1 + \left( \frac{2t}{T} \right)^2 }$$
I want to plot the frequency response of the channel as if it was a digital filter.

How can I do this in MATLAB?

I tried using symbolic variables, computing Fourier transforms of the input and the output signal using fourier command, and then finding the numerical values using subs and double commands, with no success.

I'm not sure about the symbols.
Yet if you have a function in time as the input of the channel and another function in time which is the output of the channel and assuming the channel doesn't have zeros in the spectrum and the input function as well, the channel Transfer Function is given by:

$$ \frac{{X}_{Out}(f)}{{X}_{In}(f)} $$

In MATLAB, just apply FFT on each, and do element by element division.
This will give you the DFT of the transfer function.
Apply IFFT to see the time response function.


Compute Power Delay Profile (PDP) from Delay Spread I have data pertaining to channel obtained from measurements. In order to make use of this data, I need to convert it into a form as follows:

    ChanMod.PDP_dB = [-2.5   0 -12.8  -10    -25.2  -16; % Average power [dB]
                      0 300*10^-9 8900*10^-9 12900*10^-9 17100*10^-9 20000*10^-9]; % delay (s)
    ChanMod.normH = sqrt(sum(10.^(ChanMod.PDP_dB(1,:)/10)));

However, the data from the measurements is in the form of Delay spread vs distance. I am not sure how to determine the power delay profile from the Delay spread data.

I am new to signal processing. Any help would be appreciated

Not at all. Delay Spread is just a statistical number that says something about how far the PDP is spread out. You can't reconstruct the PDP from that.

It's basically the same as wanting to reproduce the exact height distribution of a group of people when all you know is the variance of height.


Inverse Fourier Transform From Plots (2019 edition) Hello I borrowed the title for another post.

I cannot figure out how to find the inverse fourier transform from this spectrum.

I know what the transform is

I'm sorry for the plot being hard to read, it is how we were given it and I understand the frustration. The magnitude plot goes from -w0 to w0 with a amplitude 1. The phase plot has the same interval but alternates from pi/2 to -pi/2. Apologies for formatting I am new and on mobile.
Can someone help me I am totally lost.

That's actually a fairly a tricky one. Two ways you can go about it

Method 1

Magnitude is that of a ideal lowpass filter. Inverse of that is sinc function
Phase is that of an hilbert transformer (times -1). Impulse response is (roughly) $1/\pi t$
convolve the two

Method 2

Write out the equation for inverse Fourier transform
Pop in values from the graph and simplify a bit
You get either sum or integral for unity amplitude sine waves up to $\omega_{0}$

Neither one is pretty.


FFT: Does the result of $N \log_2(N)$ stand for total operations or (complex) adds and multiplies? Calculate the amount of operations you need for the Fast Fourier Transform.

Take $N=2048 \times 2048$

$N\dfrac{\log(N)}{\log(2)}$

$\dfrac{\log(2048 \times 2048)}{\log(2)}$ is 22.

$(2048 \times 2048) \dfrac{\log(2048 \times 2048)}{\log(2)}$ is 92274688.

But does $92274688$ stand for total operations or $92274688$ adds AND $92274688$ multiiplications?

So the total amount of operations is $2\times 92274688=184549376$.

{EDIT: The original question concerned a $1920\times 1080$ FFT, not a $2040\times 2048$ FFT. So the parts of this answer that deal with sizes that are not powers of two are no longer relevant. I will leave them anyway. The last part of this answer, dealing with complex vs. real operations, remains relevant.}

A 2-D FFT of size $N_1\times N_2$ is equivalent to performing $N_1$ FFTs of size $N_2$, and $N_2$ FFTs of size $N_1$. So the cost becomes
$$C_{2D}(N_1,N_2) = N_1C_{1D}(N_2)+N_2C_{1D}(N_1),$$ where $C_{1D}(N)$ is the cost of a single $N$-point, 1-D FFT. Your proposed calculations are correct, then, if you assume $C_{1D}(N)=2N\log_2 N$, because this evaluates to $$2N_1N_2(\log_2 N_1+\log_2 N_2)=2N_1N_2\log_2(N_1N_2).$$

[EDIT: This paragraph is no longer relevant; the OP has changed the FFT sizes to a power of two.] However, your other problem is that $C(N)$ is not necessarily $kN\log_2 N$ (or close to it) unless $N$ is a power of two. In particular, if $N$ is prime, then you're going to need $kN^2$ operations. For values of $N$ that are composite prime factors, then the computational complexity is a function of the number of prime factors and their values. See this page for a discussion of efficient prime factor FFTs. Roughly speaking, if the number of prime factors is sufficiently large, then the cost is going to grow like $N\log N$.

Even for powers of two, $C_{1D}(N)=2N\log_2 N$ is not quite right for several reasons. First of all, the standard radix-2 FFT actually requires $(N/2)\log_2 N$ complex multiplications and $N\log_2 N$ complex additions. Furthermore, you really don't want to just add these two counts together. After all, one complex addition requires 2 real additions, and one complex multiplication requires 4 multiplications and 2 additions. So if you want to combine these two operation counts, you need to do it in terms of these real operations, giving you $C_{1D}=5N\log_2 N$ real operations.

You can do significantly better than this, it turns out, if you take a careful accounting of "trivial" complex multiplications. One of the most efficient algorithms from a real operation standpoint is the split-radix algorithm (which requires a power of two) requires $C_{1D}(N)=4N\log_2 N-6N+8$ real multiplications and additions for a single FFT. For a 2-D FFT of size $N_1\times N_2$, the total cost would be
$$\begin{aligned}
C_{2D}(N_1,N_2)&=N_1(4N_2\log_2 N_2-6N_2+8)+N_2(4N_1\log_2 N_1-6N_1+8)\\
&=4N_1N_2\log_2(N_1N_2)-12N_1N_2+8(N_1+N_2).
\end{aligned}$$

Finally, let me add that counting floating point operations is a crude way to measure practical performance. Memory access patterns are often a more important factor in performance. Quoting Wikipedia: "In practice, actual performance on modern computers is usually dominated by factors other than the speed of arithmetic operations and the analysis is a complicated subject (see, e.g., Frigo & Johnson, 2005), but the overall improvement from O(N2) to O(N log N) remains."


Bandpass filter using wavelet transform I'm working on a speech recognition project. The first step of this project is to find phoneme in the speech signal. To do that, I found this paper that discusses about it.

In the paper, wavelets are used to visualise the signal in different frequency band. Here is my problem :

So far, I know how to decompose the speech signal using wavelet transform at different level (wavedec in MATLAB) But I don't know how to filter this signal.

With Fourier transform, a simple threshold on the FFT (focus on specific frequency band) will do the work. And as far as I understood, wavelet kinda work like Fourier so I guess it's working the same way.

If you refer to the documentation for wavedec, a signal $x$ is decomposed on one level into two sets of coefficients: $cA_1$ and $cD_1$. They correspond to a low-pass and a high-pass filter applied to $x$, following by a downsampling. As a result, if you reconstruct $x_1$ with waverec from $cA_1$ only (setting $cD_1$ to zero), and $x_2$ from $cD_1$ only (setting $cA_1$ to zero), $x_1$ will mostly correspond to the lower half of $x$ spectrum, and $x_2$ the upper half.

The same reasoning works on several levels: if your signal has a range of frequency in $[0\,,f]$, $cD_1$ gathers coefficients mostly from $[f/2\,,f]$, $cD_2$ gathers coefficients mostly from $[f/4\,,f/2]$, etc.

So for a sampling frequency of $44100$ Hertz, the bands would be:

$cD_1$: $11025 \to 22050$
$cD_2$: $5512.5\to 11025 $
$cD_3$: $2756.25 \to 5512.5$
$cD_4$: $1378.125  \to 2756.25$
$cD_5$: $689.0625 \to 1378.125 $

If you want to filter out a frequency band, you can zero wavelet coefficient whose spectrum intersect that frequency band, and reconstruct the data. This is a form of thresholding in the wavelet domain. As you guessed, it can work in a way similar  to Fourier.

Indeed, thresholding and shrinkage are very effective with wavelets, possibly more than with a Fourier transform, for denoising. In the wavelet domain, you can design the shrinkage to preserve specific time intervals, allow smooth transitions, etc.

But the wavelet filters are imperfect filters. And downsampling cause aliasing, causing a not-so-clean filtering. To perfect a pure band-pass filter, I would not recommend the DWT (discrete wavelet transform), unless the wavelet is of quite high order.


What do high and low order have a meaning in FIR filter? I'm new DSP fields, and slowly study of FIR Filer and with Matlab which is fir1 function.

My question is that

What do high and low order have a meaning in fir1 filter? I can't understand what does it have a meaning when N=2 or N=48.
In cut-off frequency, I found some explanation from the MATLAB as the below,

fir1   FIR filter design using the window method.
    B = fir1(N,Wn) designs an N'th order lowpass FIR digital filter
    and returns the filter coefficients in length N+1 vector B.
    The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0
    corresponding to half the sample rate.  The filter B is real and
    has linear phase.  The normalized gain of the filter at Wn is
    -6 dB.

B = fir1(N,Wn,'high') designs an N'th order highpass filter.
You can also use B = fir1(N,Wn,'low') to design a lowpass filter.

If Wn is a two-element vector, Wn = [W1 W2], fir1 returns an
order N bandpass filter with passband  W1 < W < W2. You can
also specify B = fir1(N,Wn,'bandpass').  If Wn = [W1 W2],
B = fir1(N,Wn,'stop') will design a bandstop filter.

But  I want to know should we have to know the sample rate to set cut-off frequency? usually in practically, How does it handle to set the cut-off frequency?

I've seen from the book, in the Liner system, if we put the impulse to the Liner system then we can totally analyze the system.
In practically, how to implement the impulse signal to the liner system to get the filter information?

1- an FIR filter can be described by the following input/output relation:
$$ y[n] = \sum_{k=0}^{N} b[k]x[n-k]$$

where $N$ is the order of the FIR filter (its length-1) and those coefficients $b[k]$ of length $N+1$ are the FIR filters impulse response, or the filter coefficients in practice.

As one can see high order FIR filter means $N$ to be large and low order means $N$ to be small. Note that most of those coefficients can be set $0$ but the last one at $x[n-N]$ would still make a high order filter (specifically a filter with a pure delay).

The delay associated with FIR filters (of symmetric impulse response types) is directly related to the order $N$ of the filter.

2-
An FIR low-pass filter is a mathematical algorithm and its cutoff frequency is normalized to be between $0$ and $\pi$ radians mathematically and between $0$ and $1$ for software such as MATLAB. When you want to find the physical cutoff frequency of such an FIR filter , you should consider the physical sampling rate and relate it through $$f_c = \frac{\omega_c F_s}{ 2\pi} $$
where $\omega_c$ is the cutoff frequency (in radians per sample) of the mathematical filter, and $f_c$ is the physical cutoff frequency in Hertz.

3- Impulse response of the FIR filter is directly available in the coefficients $b[k]$ as I've explained. When you are dealing with an IIR filter or an LTI filter which does not posses a LCCDE expression, then practical computation of the impulse response that can be done in a number of ways. The simplest could be to create a signal that defines the unit impulse as delta = [1, zeros(1,L-1)] then the $L$ points of the impulse response can be computed directly by h = filter(b,a,delta) . Note that if impulse response $h[n]$ does not decay fast enough then you should increase the length $L$.


Music source separation: Incomplete notation in published IS-NMF/EM pseudo-algorithm? In a paper by Fevotte et al. on IS-NMF, which describes algorithms/implementation for multiplicative and EM variants of IS-NMF, the following algorithm is provided for IS-NMF/EM:

In the algorithm's footnote, the authors state that the VHAT matrix $\hat{\mathbf{V}}$ only needs to be calculated during initialization, but then go on to say that $\hat{\mathbf{V}}$ is subsequently updated as $\hat{\mathbf{V}}=\mathbf{W}_{-k}\mathbf{H}_{-k} + \mathbf{w}^{new}_k h^{new}_k$ during each iteration.

Therefore, I believe that the hat matrix $\hat{\mathbf{V}}$ actually needs to be updated inside the iterations.  So my question is: Before I start coding anything, shouldn't an update for $\hat{\mathbf{V}}$ be appropriately placed inside the iterative loops?  Also, where would you place it?  Better yet, where would you place all the updates mentioned in the footnote, so that the footnote is not needed.  (why mention updating of coefficients/parameters in a footnote if they can be placed in the algorithm?)

The implementation in MATLAB is here:
https://www.irit.fr/~Cedric.Fevotte/extras/neco09/code.zip


Frequency Components of an FFT I was reading a couple of these threads and I was wondering:
Is it possible for an FFT signal to include frequencies less than $1/N$?

Consider a time-series of length $N$, and I perform an FFT on it. In the zero-th bin of the FFT is the DC signal, and first bin is frequency $1/N$, in the second bin is frequency $2/N$ etc.

Does this mean that the smallest frequency accounted for is a signal where one full period fits in the data length? Or can there also exist frequency where only half/quarter periods will fit in the given length?

A DFT will convert a segment of a periodic waveform of any frequency (uniquely below Nyquist)  into an orthogonal set of frequency results.  The problem is that if the signal is not periodic in the FFT length, the result will not be in one result bin, but be smeared into all the FFT result bins, due to convolution with the transform of the rectangular or other window function, making it difficult to separate the periodic signal from narrowband + wideband noise (or DC + ramp + noise for spectrum below the 1st bin).

Another potential problem with estimating non-periodic real spectrum very near DC (or Fs/2) is interference from its own negative frequency image in the DFT/FFT result, which can be constructive or destructive depending on phase, and thus non-trivial to account for.

So yes, the low frequency, below Fs/N, will "exist" in the DFT/FFT complex result, but be more difficult to impossible to find or estimate, depending on the signal to noise ratio.  But in zero noise, you theoretically need only 3 or 4 non-aliased samples to find a single pure sinusoid of any frequency, so it will be represented in the transformed N samples, for N >= 4, just not clearly in only 1 or 2 DFT/FFT result bins.


Correlated fading channel: transmit and receive antenna correlation matrices I am trying to create a non-i.i.d channel: correlated Rayleigh fading channel.
$$ H= R_H \cdot H_o,$$
where $H_o$ is i.i.d. channel ( uncorrelated) and $ R_H$ is a channel correlation matrix, is defined as
$$ R_H= R_T^T \otimes R_R,$$
with $R_R$ receive antenna and  $ R_T$ transmit antenna correlation matrices.

My question is how to find $R_T$ and $R_R$, if numbers transmit and receive antennas are given only?
Can I assume them as identity matrices?

$R_H$ could be anything other than identity and you would get correlated fading with the definitions you have in place in the question. Defining a kronecker matrix product is one way.

$R_H$ should be $N_R \times N_R$ matrix, you can make the kronecker product to result in a non identity matrix of this size.


Can we normalize data to unit variance without distorting the phase? Consider we have one dimension time samples. The noise in the samples made some portions of the time samples amplitude very high. Obviously the short variance for those time samples will have high variance. I want to do running variance (of some window size)  and normalize data by dividing by the variance. Is this procedure would destroy the phase of the phase of the signal? Is there a way to normalize data with unit variance without distorting the phase?

"Normalizing the variance": Since variance is a real number, normalizing it means multiplying the signal with a real number.

That doesn't change the phase, ever.


What is the average of $\DeclareMathOperator{\rect}{rect} \rect(\cos(\pi t/2)) $? We have this signal:
$$\operatorname{rect}\left(\cos\left(\frac{\pi t} {2}\right)\right) $$
I must find the average power , how can i get there ?
My solution:
I have seen that
$$-\frac 12 < \operatorname{rect}\left(\cos\left(\frac{\pi t}{2}\right)\right) < 1$$
is for
$$\frac 23+k \pi <t < \frac 43 + k\pi$$
so in a period it should be rectangular pulse , but how do I calculate the average ?

The given signal is clearly periodic. So just follow these steps and find the solution:

Find the period $T$.
Figure out the interval within a period for which the signal equals $1$.
Compute the average: $$\overline{x(t)}=\frac{1}{T}\int_0^Tx(t)dt$$

If you're convinced that you did everything right, don't worry about a given solution which is different.


Clock Jitter and Phase Noise I would like to know the impact of the time jitter of the clock driving an ADC on the appearing phase noise while sampling a sine wave. Of course, with strong the clock jitter, the sampled sine wave seems to exhibit more and more phase noise. But how much ?
  Regards,
     Mike

For a sine wave with a period $t_p$ the clock jitter or clock edge time error $t_j$ gives a phase error of $2 \pi t_j/t_p$ radians, where $2\pi$ represents a full period of a sine wave in terms of phase. The period $t_p$ of the sinusoid is an inverse of its frequency $f$, so the phase error can also be written as $2 \pi t_j f$. No matter if error is the peak error or the root mean square error, or whatever, as long as it is the same kind for both clock edge time error (jitter) and phase.


How to classify the movement type of a target into two categories(direct motion or else) by using a doppler radar? How to classify the movement type of a target just into two categories(one is direct motion which is almost straight line movement and the other is anything else) by using a doppler radar with a microdoppler technic, or the spectrogram of the doppler frequency?

Could you recommend a paper, patent, and something like that related to this issue?

Thank you so much.

If range isn’t important, sure you could do something like this with unmodulated pulses (i.e the pulsed you transmit are just sinusoids carriers with some known frequency). You would transmit the pulse and then receive it at some time delay. You could attempt to construct some sort of spectrogram or other time-frequency representation of the received signal. If the target is moving at a constant velocity, there would be a constant Doppler shift. If the target had non-constant velocity, you would be able to observe a non-constant Doppler shift.

If you cared about resolving range to the target, you would have to use a pulsed Doppler system. That would be very difficult to do with a pulsed Doppler system. Pulsed Doppler radars tend to ideally assume the target signal is moving at a constant rate, and that over the entire coherent processing interval (CPI), it doesn’t change range bins. The reason the assumption is key is because if the target shifts range bins, you will not get all of your Doppler processing gain and you’ll tend to smear out the target on the range-Doppler map, which makes detection of the target difficult if not impossible

In industry, the issue of range cell migration over a CPI for high velocity targets is pretty difficult but doable under the right circumstances. You can check out something called the keystone algorithm, that has seen some use in practice. I’d say it’d be very important to select your bandwidth and sample rate such that you don’t cause any extra apparent range migration of the target.

As for microdoppler on pulse Doppler radars , I have yet to actually see that be used in the industry I work in (defense), and that community largely treats it as computationally inefficient. With that being said, let’s say you have all the computation time in the world and you could do such a thing, sure it may reveal some characteristics. One classic example would be being able to detect the type of engine on some sort of airplane based on the microdoppler characteristics, e.g. the rotors of a helicopter or the supersonic wake of a jet.

In summary: yes, if you didn’t need range information you could analyze Doppler through using a spectrogram. Keep in mind that your RF carrier frequency will determine the Doppler shift that is apparent, so selecting an appropriate RF is pretty key for something like this.


Extracted envelope via Hilbert Transform is tracking the signal I am trying to extract the envelope of the signal shown in blue, which is the low-pass filtered output of a second-order system. Low-pass filtering was done in order to reduce measurement noise. To be able to determine the decay of the system I want to extract the envelope of the signal. My MATLAB code looks as follows:

d = fdesign.hilbert(500,0.05);
hd = design(d,'equiripple','SystemObject',true);
coeff = hd.Numerator;
h = coeff / sum(abs(coeff));

y_hilb =  1i * filter(h, 1, y_filt);
y_an = y_filt + y_hilb;
env = abs(y_an);

Unfortunately, the envelope that I am extracting is actually tracking the signal rather than showing the actual envelope. What am I doing wrong here?

Sampling Period: 0.01s

EDIT1:
I changed the code as follows. Unfortunately, the result looks similar.

d = fdesign.hilbert(501,0.05);
delay = (N - 1) / 2;

y_an = y_filt(1 : end - delay) + 1j * y_hilb(delay + 1 : end);

EDIT2:
As suggested, I removed the scaling of the coefficients. Result:

Computing the envelope of a signal via the Hilbert transform only works well for bandpass signals.

But apart from that, you're not taking into account the delay of the Hilbert transformer when creating the analytic signal. It would be easier to choose an odd filter length $N$, because then you get an integer delay:

$$d=\frac{N-1}{2}$$

Also note that there's the command hilbert for directly computing the analytic signal. Furthermore, there's a function envelope that lets you choose different methods for computing the envelope of a signal.


Sinc interpolation of pure sine wave sampled just at Nyquist frequency Following this question: Shannon-Nyquist theorem reconstruct 1Hz sine wave from 2 samples

could you explain the algorithm to apply for sinc interpolation to avoid the "sawtooth" effect due to linear interpolation?

(It seems to me that the https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula  Shannon-Whittaker formula would be the one suitable for this?)

Here again let me note that exact Nyquist frequency for a pure sine wave should be avoided. Shannon-Nyquist sampling theorem requires that there's no impulse at the Nyquist frequency as a consequence of bandlimitedness, the content at the exact Nyquist frequency is taken to be zero.

Then the following code demonstrates the approximate simulation of an ideal sinc based interpolator applide to (near) critical samples of a pure sine wave. Note that, any finite observation of a signal cannot be bandlimited, so this simulation is not a perfect representation of the true output from an ideal interpolator, nevertheless by chosing signal duration long enough, one can attain an approximately bandlimited signal.

f  = 1;                % 1 Hz. sine wave...
Fs = 4.2*f;            % sampling frequency Fs = 2.2*f  ; a bit more than the Nyquist rate.
Td = 25;               % duration of observation ultimately determines the spectral resolution.
t  = 0:1/Fs:Td;        % observe 25 seconds of this sine wave at Ts = 1/Fs
Td = t(end);           % get the resulting final duration
L  = length(t);        % number of samples in the sequence
M = 2^nextpow2(10*L);  % DFT / FFT length (for smoother spectral display, not better resolution! )

x = sin(2*pi*f*t);     % sinusoidal signal in [0,Td]
%x = x.*hamming(L)';   % hamming window applied for improved spectral display

% Part-II : Approximate a sinc() interpolator :
% ---------------------------------------------
K = 25;                  % expansion factor
xe = zeros(1,K*L);       % expanded signal
xe(1:K:end) = x;

D = 1024*8;
b = K*fir1(D,1/K);     % ideal lowpass filter for interpolation

y = conv(xe,b);
yi = y(D/2+1:D/2+K*L);

subplot(3,1,1);
plot(t,x);
title(['1 Hz sine wave sampled at Fs = ',num2str(Fs),' Hz, Duration : ', num2str(Td), ' s'])
%xlabel(' time [s]');

subplot(3,1,2);
plot(linspace(-Fs/2,Fs/2-Fs/M,M),fftshift(abs(fft(x,M))));
title(['magnitude of ', num2str(M), '-point DFT / FFT of y[n]']);
%xlabel('Frequency [Hz]');

subplot(3,1,3)
plot(linspace(0,Td,length(yi)),yi);
xlabel('approx simulation of ideal sinc interpolation');

With the result of


How to detect a signal while ignoring another signal which is a phase shifted version? I have a signal coming from a sensor as shown below. This signal has high frequency noise (as can be seen from small fast wiggles) and four regions of interest where distinct spikes can be seen. The first region of interest starts with signal dropping first and settling. Second interesting region shows similar pattern; it drops first. The third region of interest is different from first and second because this time it spikes upwards and then settles. And fourth area of interest is also similar to third one; it spikes upwards first and then settles.

I consider third and fourth area of interests as similar signals as first and second but phase shifted versions.

My goal is not only to detect these interesting areas but also distinguish between the first two and last two. I then designed a matched filter and passed the raw signal through this matched-filter. The template for the matched filter was designed to detect the first two variants (the ones which fall downwards first) as shown in the picture below. Note that template is time-reversed version of the first area of interest in the signal.

The output from matched filter is shown in the figure below which improves the SNR quite a bit and allows for easy detection of areas of interest. As can be seen, I can very clearly see the spikes corresponding to the four areas of interest. However, I cannot distinguish between the first two and second two. How do I do that?

Try this: Apply a moving-average-window to the output signal. A window with a length close or equal to your matched-filter template signal. This window is supposed to add all the samples within its range. When this window applies to the location of first and second waves, the result's sign would be different to that of third and forth. You can also try a weighted window like Hanning window.


Python code for QPSK modulator Have any of you done anything like QPSK Modulation in python? I am in the middle of it and having problems in setting up bits? like odd and even bits.

Have you got any idea?

Regards,

Sam

Are you using scipy by any chance ?

If yes , this might help

https://scikit-dsp-comm.readthedocs.io/en/latest/_modules/sk_dsp_comm/digitalcom.html

check for functions "QPSK_tx"  and "QPSK_bb"


Help required with BM3D Image denoising algorithm implementation This is what I am doing and will be grateful for any help.

I selected 16 similar blocks to a reference block from a search block, put them in a 3D array, Wavelet Thresholded them, took inverse transform.
These 16 similar blocks belong to one 32 * 32 search block as mentioned above.
Now I should multiply each of these blocks with their weight, and then add them all up, creating a 2D array from a 3D array.
This new block then replaces all the 16 original blocks that I had selected before for hard thresholding.

Am I correct here? I will be absolutely grateful for a response.

Step 3 is not correct. The first part is ok (each denoised result is weighted by some likelihood of correctness, that was obtained previously from the variance of the data), but you don't have to add them up. Each slice from the 3D array is a 2D patch that is back-projected to the resulting image.

The sum happens when 2 back-projected patches do overlap in the final image domain.

About back-projection

Back-projection is the process of taking a 2D slice from the denoised stack and put it back in image space.
This is required because unlike most variants of NL-means BM3D denoises patches and not the central pixel of a patch.

So, you have 2 images when implementing BM3D:

the original (noisy) one;
the result (denoised) image.

The spatial location of a patch should not change during the denoising process, otherwise the image content would be destroyed. So, each slice from the denoised patch stack goes back to its original place.

Now, one issue remains: what to do when back-projected patches overlap? What value should we keep for a single pixel?
The solution chosen by BM3D is that each pixel contains a weighted mean of the patches that overlap in each pix.

To implement this, you can simply pre-multiply each slice of the denoised stack by its weight, add a third image to store the sum of the weights inside each pixel, sum for each pixel its value from all the denoised patches that contain it, and in a final step (once all the patch stacks haves been processed) divide the result image by the weight image.

Image processing main loop (outline)

There is no "unleft" pixel in the final result. The outline of the main loop is as follows:

for every pixel of the image, form a reference patch;
for a given reference patch, find n <= N (B=16 in your case) patches whose distance to the reference patch is less than some threshold and make a stack;
denoise the stack;
back project the patches in the denoised stack;
back to 1) for the next pixel).

Thus, you always have one stack for a given patch (degenerate case: only 1 patch in the stack), and a pixel is always involved in several patches.


Why do we use white noise to estimate echoes? I'm interested in the case where I have received a signal as a composition of three echoes of the transmitted signal $x[n]$ (each echo has its amplitude and time delay). I know that in order to compare the signals I use cross-correlation. My question is if I use $r_{xx}[m] = \delta[m]$, would that make it easier to calculate the parameters of echoes? And, if so, can you explain me why?

The mathematical justification would be that the Cramer Rao Bound on the estimate of time delay is inversely proportional to the time bandwidth product of the pulse, so in cases where there are many distinct echoes, the shorter the pulse is better to resolve each arrival, so one would increase the bandwidth to minimize the bound on time delay.

In the case where the signal is long, as illustrated by the matlab code below for a single delay:

  clear all
  x=randn(1,16*32768);
  x1=[x zeros(1,.5*128)];
  x2=[zeros(1,.5*128) x];
  figure(1)
  pwelch(x1+x2)

The psd has periodic nulls, and you can get an estimate from the null spacing, and the wider the band (white noise), the more nulls that can be resolved.

$$
H(\omega)=1+\alpha e^{\jmath \omega \tau}
$$

For 3 delays it will be a bit more of a bother, but doable.  Actually if you don't know the waveform, like a passive SONAR problem, this is an essential approach.