How to make spectrogram smooth? I am quite new in DSP and generating the spectrogram of an audio file. My spectrogram is not smooth and it is showing the quite raw image with pixel values, something like this

While I am looking for a smooth spectrogram like this

Where I am doing mistake? Is it because of size of window size? My code to generate Mel spectrogram is
def readData(file):
    origData,origSampFreq = librosa.load(file, sr=None)
    return origData, origSampFreq

def resample(originalData, origSampFreq, targetSampFreq):
    resampledData = librosa.resample(originalData, origSampFreq, targetSampFreq)
    return resampledData

def normalizeSound(resampledData, axis):
    """ Axis is 0 for row-wise and 1
    for column wise"""
    normalizedData = normalize(resampledData, axis)
    return normalizedData

def calculateMelSpectogram(normalizedData, hop_length, win_length, sr):
    #newSamplingFreq = 16000
    S=librosa.feature.melspectrogram(normalizedData, sr=sr, hop_length=hop_length, win_length=win_length)
    return S

#Plot melspectogram

def plotMelSpectogram(S, sr, name, ref=np.max):
    plt.figure(figsize=(10,3))
    S_dB = librosa.power_to_db(S, ref=np.max)
    librosa.display.specshow(S_dB, x_axis='time',y_axis='mel', sr=16000,)
    plt.colorbar(format='%+2.0f dB')
    plt.title('Mel-frequency spectrogram')
    plt.savefig('./chunk_images/' + name + "mel.png",dpi=(300), bbox_inches='tight')
    plt.tight_layout()
    plt.show()
def featureExtraction(audioFile, name, targetSampFreq = 16000,
                      axis =0 ,
                      hop_length= 256,
                      win_length=512):

    y, y_sr = readData(file=audioFile)
    print(y, y_sr)
    resampledData = resample(originalData=y, origSampFreq=y_sr, targetSampFreq=targetSampFreq)
    normalizedData = normalizeSound(resampledData, axis=axis)
    S = calculateMelSpectogram(normalizedData=normalizedData, hop_length=hop_length, win_length=win_length, sr=targetSampFreq)
    plotSound(soundData=normalizedData, sr=targetSampFreq,x_axis_string='time' , name = name)
    plotMelSpectogram(S, sr=targetSampFreq, name = name, ref=np.max)
    return S

# plot orginal time domain data

def plotSound(soundData, sr, x_axis_string, name):
    plt.figure(figsize=(10,3))
    waveplot(soundData, sr, x_axis=x_axis_string)
    plt.savefig('./chunk_images/' + name + "sound.png",dpi=(300), bbox_inches='tight')

Try increasing the overlap in the spectrogram calculation, and using an interpolated array-to-image renderer (shading interp in matlabesque).
-k


What are the advantages of Laplace Transform vs Fourier Transform in signal theory? What are the advantages of Laplace Transform vs Fourier Transform in signal theory?

Laplace transforms can capture the transient behaviors of systems.  Fourier transforms only capture the steady state behavior.  Of course, Laplace transforms also require you to think in complex frequency spaces, which can be a bit awkward, and operate using algebraic formula rather than simply numbers.

If you want to see the power distribution of a signal over time, Fourier transformations are often the easiest way to do it.  However, if you want to understand what a system does when you flip a light switch, you typically need Laplace transforms.


Save Matlab audio output I need to save the audio output of a script in Matlab but I don't know how to do it. It is a complex sound, with several pauses of silence with different durations, so I cannot simply use audiowrite (for my understanding, I need one single matrix S to save a file using audiowrite(S,Fs)).
A simple example of how the sound is generated would be the following:

for ii = 1:length(cf) - 1
sound(s(ii,:), sf);           % sound presentation
pause(d + 0.025);             % waiting for sound end + silent period
end

where d is the sound duration and 0.025 is the duration of the silence.
I was thinking in creating a vector describing the silent period (similarly to the matrix s describing the sound), but as it has a different duration with respect to the sound so the number of columns (which represents the number of signal samples) would be different.

Is it possible to record the output? I would like to save exactly what I hear in the speakers after running the script. It is my first time with sounds in Matlab and I did not find anything online. Thank you for your help.

Are you aware of cell-arrays in Matlab?

Consider the following short script-file,

    close all; clear; clc;

    x=randn(1024, 1);
    y=2+sqrt(0.5).*randn(2048, 1);

    S=[ {x}; {y} ]; size(S)

    z=[ S{1, :}; S{2, :} ];

    figure; ...
        subplot(3,1,1); ...
            plot(S{1, :}); grid on; xlim([1 numel(S{1, :})]);
        subplot(3,1,2); ...
            plot(S{2, :}, 'r'); grid on; xlim([1 numel(S{2, :})]);
        subplot(3,1,3); ...
            plot(z, 'k'); grid on; xlim([1 numel(z)]);
        shg;

Once you create a vector of cell-arrays (each cell-array will have a different length that is dependent on the duration of the silent periods), you can save it to a file.

I hope this helps.


Plotting a Time Frequency contour/colormesh plot of a Discrete Wavelet Transform I have a pressure vs time data of a noise on which I wish to perform discrete wavelet transform.
I have divided my frequency range into 1/3rd Octave Bands and have calculated sound pressure level at each band.
I am very much confused on how to perform a discrete wavelet transform in 1/3rd Octave Frequency Bands or any set of bands for that matter and further how to interpret the output(coefficients) so as to plot a Time-Frequency colormesh or a contour plot of my results.
I was able to plot a continuous wavelet transform of my data using the pcolormesh module of matplotlib. The output of pywt.cwt are easy to interpret and hence I could plot a colormesh for it.
I am a beginner in DSP and have just started learning wavelet transforms, so any help will be much appreciated. Thankyou.

With wavelets, that rely on scaling, you can expect time-scale, more than time-frequency. Obtained 2D images of scalograms from continuous wavelets is relatively easy. Most of them are not exactly invertible.
Exactly invertible discretized wavelets come in different shapes. To preserVe images with equal-sized pixels, you could prefer redundant discrete wavelets over DWT.
So, maybe rwt.rdwt()


Need help with FFT results I am playing around with FFT and cant seem to get the expected results. I have a sound clip that i recorded from a mic placed on an analogue clock. So it has a sound every second.

[wave,fs]=wavread('clock.wav');
t=0:1/fs:(length(wave)-1)/fs;
plot(t,wave);
n=length(wave)-1;
f=0:fs/n:fs;
wavefft=abs(fft(wave));
plot(f,wavefft);

This is the code that i am using. And here are the results.

I should be seeing a huge spike at 1hz but i am not.

The FFT followed by the sound wave.

Since you are looking for a low frequency periodicity in the amplitude envelope of your sound clip, and not the tone waveform or pitch, you might try squaring or taking the absolute value of your signal before applying the fft, and then ignoring higher frequency bins in the result (similar to a cheap AM demodulator).


Multiple-image dense point cloud reconstruction with camera extrinsics/intrinsics I have the following data available:

Images from different angles
X, Y, Z, pitch, yaw, and roll for each camera/image, together with the timestamp
Lens distortion parameters and intrinsic camera parameters such as focal length and CCD width

On my dataset, OpenSfM produced poor reconstructions and I realized that it tries to figure out camera extrinsics, which I already know.

My question: I am trying to produce a dense point cloud. Is there a name for this problem when I have such prior information? Which publicly available libraries can help me accomplish this task? Can someone point me in the right direction?

Thanks!

First, a warm welcome to SE!
Basically, you have a calibrated 3D reconstruction problem. The typical approach follows a 5-stage pipeline:

Identify 2D features in each image along with the associated descriptors. Algorithms such as SURF, SIFT or AKAZE are heavily used and are available in many vision libraries such as OpenCV.
Match the extracted keypoints across the whole dataset by using their descriptors. Once again, this is a well studied topic. The keypoints and their multiview association will be called feature matches and will form the basis for 3D reconstruction.
Extend the pairwise matches to the multiview case by linking the cameras to one another, associating the ones that have enough common points. In MATLAB, the relevant function would be findTracks.
Use an N-view triangulation algorithm (typically a multi-ray intersection) to find an initial estimate of the 3D points. In MATLAB one could use the triangulateMultiview to solve for the multiview intersection. Normally, prior to this stage, one would solve for the camera poses, but since you have it already, I'm skipping this part.
Now use bundle adjustment to refine the 3D structure. I am not sure how you obtain the poses, but it is strongly recommended to also adjust the camera extrinsics additionally to the points. Libraries such as Ceres Solver are capable of solving for both structure and motion, simultaneously.

Optionally, a post-processing step consisting of densification can be run. PMVS could be of help here:

  Accurate, dense, and robust multi-view stereopsis. Yasutaka Furukawa
  and Jean Ponce. IEEE Trans. on Pattern Analysis and Machine
  Intelligence, 32(8):1362-1376, 2010.

It is another possibility to benefit from existing multi-stage SfM pipelines, such as OpenMVG. You would simply replace the OpenMVG poses by your camera poses and keep the rest of the workflow the same. This way, you get some outlier treatment and post-processing steps for free. MATLAB also has an SfM toolbox that you could easily test.


What is the right way to calculate the inverse Z-transform of $zX(z^{-1})$ say the signal $x(n)$ has the z transform $X(z)$ and there is signal $x_1(n)$ that

$X_1(z)=zX(z^{-1})$

I tried 2 different approach to get the relationship between $x(n)$ and $x_1(n)$ and the result is different without any error being noticed

Approach 1:

$x_i(n)=x(-n) \rightarrow X_i(z) = X(z^{-1})$

$x_1(n)=x_i(n+1) = x(-n-1) \rightarrow X_1(z) = zX_i(z)=zX(z^{-1})$

Approach 2:

$x_i(n)=x(n-1) \rightarrow X_i(z) = z^{-1}X(z)$

$x_1(n)=x_i(-n)=x(1-n) \rightarrow X_1(z) = X_i(z^{-1}) = zX(z^{-1})$

Hence in the end, I don't know $x_1(n) = x(-n-1)$ or $x_1(n) =x(1-n)$ and I can't find any mistake. I know that different function can have same z-transform but different ROC but in this case, no any info of ROC is given

The first result is correct, the second is wrong. You're right that $x_i[n]=x[n-1]$ corresponds to $z^{-1}X(z)$. And if you replace $z$ by $1/z$ you need to replace $n$ by $-n$. But $x_i[-n]$ is simply $x[-n-1]$,  and not $x[-n+1]$. Note that you don't replace the whole argument of $x[n-1]$ by its negative value, but you just replace $n$ by $-n$.


Affine transformation without imwarp I program function to calculate the overlap area between 2 images (which is called fixed and floating) after transformation but it takes much time and it's complex for me to combine more than one transformation parameter together.

In Matlab, the built in functions affine2d can make affine transformation (rotation, translation, and scaling) that can solve my problem. However, this function depends on imwarp which fills the empty space with new pixel (default: black).

As you can see the below image, my ROI inside the red box in the both images.

How can extract this area from the transformed image (float)?

Update: the hardest part is to find the ROI coordinator in the transformed image.

To find coordinate or index of drew region, first you have to find the coordinate of 4 colored dots, which is easily obtained by using the transformPointsForward function over the intrinsic coordinate of these 4 points in your original image. intrinsic coordinate is the index of your image matrix.

then you have to do some geometric calculations. first, find the length of one of red lines which are the (length of full warped image)/(sin(theta) + cos(theta)), then distance of red line to the nearest colored point is (length of red line)*cos(theta)*sin(theta). theta is the rotation angle. using this distance and the coordinate of 4 colored points you could find the coordinate of your ROI corners.


Noise Reduction Moving Average FIR I read that the noise reduction of a moving average FIR can be expressed as the square root of N where N is the length of the filter.
How do you prove? I can't quite get it.

Since this is homework/a learning assignment, I think a hint is worth more than the (admittedly done on this site a couple of times) complete proof:

You forgot to also read the context: That's not true for just any noise and any signal, it's only true for noise that's white and a signal that fits perfectly into the filter shape.

Think about what "white" has to do with a moving average. Hope this helps.


Frequency domain to time domain conversion-(varying sample count) I need to convert 150 Frequency domain samples to the time domain. The problem is, I have been provided with 250 time domain samples for comparing my results.  How am I supposed to relate the 150 samples to 250 samples? I can't find any relationship like interpolation/decimation between the two.

Does this have anything to do with windowing (hamming,hanning,etc), sampling frequency, min and max frequency in the original signal?

--

@Jim: I have 150 frequency domain samples and i need to get 250 time domain samples using them. If I do a normal iFFT, I would get only 150 TD samples. So I tried to append 100 zeros to the FD samples and then did the iFFT. But then, my samples are different from the original answers.

The input set I have is : 150 TD samples, sampling frequency=10GHZ, Min and Max frequency of the original signal = 500MHZ and 3050 MHZ, Window= Kaiser (beta=6)

I believe that I have been provided with either insufficient or incorrect data.
Am i right?

Consider the DTFT of an N-sequence $x_n$:

$$
\begin{align*}
X(e^{j\omega}) &= \sum_{n=0}^{N-1} x_n e^{-j\omega n} & \text{(DTFT)} \\
&= \sum_n \left ( \frac{1}{N} \sum_{k=0}^{N-1} X_k W_N^{-kn} \right ) e^{-j \omega n} & \text{(subs. IDFT)}
\end{align*}
$$

Here $X_k = X(e^{j\frac{2\pi}{N}k}), \, W_N = e^{-j\frac{2\pi}{N}}$, and let $\omega_k = \frac{2\pi}{N} k$, then rearrange:

$$
\begin{align*}
X(e^{j\omega}) = \sum_{k=0}^{N-1} X_k \left ( \frac{1}{N} \sum_{n=0}^{N-1} e^{j(\omega_k - \omega) n} \right )
\end{align*}
$$

Hence the DTFT of an N-sequence is uniquely determined by N frequency samples uniformly spaced between $0$ and $2\pi$. The term in the parentheses is an interpolation function (it's actually the digital sinc).

Check whether your $x_n$ is periodic to 150 samples.

What happens when you zero pad the DFT to length $P > N$ (add $P-N$ more $X_k = 0$ terms), and take the IDFT? You interpolate to a higher sampling rate (since $\omega_k = \frac{2\pi}{P}k$ is less than before).

Note that you're adding points after the first $N/2$ samples, explained here.


draw sines at higher freq I'm trying to plot some sine waves (code example in plain js here).
When freq is "low" (freq = 10 hz in that case), the plot is quite nice:

The problem is when I increase the freq (try to set var freq = 50 for example):

lots of ripples, it becomes distorted and not so good as plot. If I increment it more, even worse (var freq = 8030 for example is terrible).
When I see those kind of graph on pro systems, they are displayed just fine.
How would you improve it? FFT, splines, whatever? Which is the right approch?
I don't really need accurancy (i.e. for waveform analysis or whatever), just plot it nicely (as in Desmos https://www.desmos.com/calculator/eodkjlywjh, for example).

For a given «visual accuracy», you need to sample the sine at a sufficient number of time-steps per period.
At some point the display pixel density will be to low to render a sine accurately.
For waveform editors (that usually have access to only a finite rate sampled waveform, not a continously defined trigonometric function),  I assume that some practical approximation to the sinc(x) function is used to interpolate those samples to the desired density. There is still a question of what to do when the display is the limiting factor, rather than the input.


How do I get excatly the same peak amplitude versus frequency when performing FFT? I was working on simulation of steam turbine vibration under unbalance condition.
I wish to perform an FFT for the vibration velocity (mm/s) respond and plot the frequency on the x-axis and the real amplitude  on the y-axis. The respond vibration velocity data can be downloaded here:
https://drive.google.com/open?id=0B21C655TFYWJSVFjY1VPeVg0WHM

I am using the following MATLAB code:

Fs =1000; % Sample Frequency
ti = 0:0.001:15;

ti= ti';
data = signal;
N = length(data);

Y(1)=0;
Y=fft(signal);
freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = abs(Y(1:floor(N/2)))/floor(N/2);

subplot(211)
plot(ti,data)
xlabel('Time[s]')
ylabel('Amplitude')
title ('Time domain signal')
grid on

subplot (212)
plot(freq,amplitude);
xlabel('Frequency [Hz]');
ylabel ('Amplitude');
title('Amplitude Spectrum')
grid on

when executed the code this is what I've got

The time in "time domain signal" was corresponding to increased the rotation load of the rotor . For example at time 0-2 second the rotor RPM was increasing , so the unbalance force will excite the system which is lead to greater amplitude. after the 2s the rotor has reached the maximum load , so the amplitude decreased because the presented of damping inside the system
from the time domain the max peak was in 39.42 mm/s at time 1.248, and then the amplitude goes down and stable with average in 14 mm/s at the rest of the time. I Thought it will give me two spikes at frequency domain with exactly the same amplitude in y-axis which is given above

Why I didn't get the spectrum signal as I wanted, Can someone please help me? Thanks.

The reason why you do not get "...two spikes at frequency domain with exactly the same amplitude in y-axis which is given above...", is because the Discrete Fourier Transform attempts to decompose the given signal as a whole.

Given the explanation about the origin of the signal, it is very similar to that of a chirp and the chirp's spectrum is not too far away from what you get, provided of course that one accounts for the envelope shape and the fact that in your waveform, the steady state is much longer than the transient (spooling up) and this explains the strong spike at 50Hz.

If you wanted to track the amplitude of the oscillation that is caused by the imbalance with respect to time (or load), then you would have to use the Short Time Discrete Fourier Transform which is nothing more than evaluating a shorter length Discrete Fourier Transform over (potentially overlapping) windows. For example, given some $x[n]$ with $N$ samples, you could $\mathcal{F}(x[0 \ldots 255]), \mathcal{F}(x[128 \ldots284]), \ldots $. This would run a series of 256 point FFTs with 50% overlap.

Eventually, this brings us to the spectrogram. The spectrogram is the depiction of the above process where the $x$ axis now corresponds to time and the $y$ axis corresponds to frequency. The "colour dimension" represents the amplitude of a given frequency present in a signal frame.

So, in the case of your signal and assuming a constant amplitude "disturbance" across the full range of RPM, you would observe a Spectrogram with a ramp in the beginning (spool up) followed by a constant line. The slope of the ramp is equal to the rate by which the rotor spools up.

However, the "disturbance" here is not constant and it presents resonance at specific frequencies / RPMs. Therefore, if you were to run a spectrogram on this particular signal here, you would still get the ramp but with a very intense low end (at the resonant frequencies) followed by a constant amplitude line at the steady state.

For more information on how to run a Spectrogram in MATLAB, please see this link.

Hope this helps.


Calculating the Norm of Each Sample of the DFT Trasform Given a discrete signal $ \left\{ x[n] \right\}_{n=0}^{N-1} $ and its Discrete Fourier Transform $ X[k] = DFT \left\{ x[n] \right\}[k] $.

I'd like to effectively compute the norm of each sample of the DFT, namely  $ \left | X[k] \right | ^ {2} ,  k \in \left\{ 0, N-1 \right\}$ as opposed to the norm of the DFT series (Which equals the norm of the original series).

Could that be calculated without using complex multiplications (Since it's real)?
Namely, I'd like to not have to go through the calculation of the DFT and then calculate the norm of each sample.
What if $ x[n] $ is also a symmetric signal (Auto correlation for that matter), I assume then using the Cosine part makes it possible?

Thank You.

The simple answer: yes this can be done in various different ways:

Calculate the $a_{k}$ and $b_{k}$ coefficients of the regular
Fourier Series (http://en.wikipedia.org/wiki/Fourier_series) and then $a_{k}^{2}+b_{k}^{2}$ is the norm at frequency k. This is equivalent of splitting the original DFT equation into it's real and imaginary parts
Calculate the auto-correlation of $c_{xx}[n]$ of $x[n]$. Then calculate only the cosine coefficients $a_{k}$ of the Fourier Series of $c_{xx}[n]$. $a_{k}$ is the norm at frequency k.

However, each of this methods is computationally much less efficient than simply doing the DFT and then taking the magnitude squared. Using the DFT (with complex math) is by far the most efficient way of doing this. If you worried about computationally efficiency, you can shave about 1/3 to 1/2 of the computation time by using the fact that your input signal is real. The DFT of a real function can be calculated by using DFT of a complex signal of half the length. That requires a bit if extra math, which I'll skip here.


If the frequency response of a system is time-independent, is this system LTI? Intuitively, I feel like a time-variant system would necessarily have a time-dependent frequency response, and vice-versa. So, is the time-independence of the frequency response necessary and sufficient to say that a sytem is LTI ?

LTI means Linear Time-Invariant systems. The system has to satisfy two conditions.
(1) Linear and (2) Time-Invariant.
(1) Linear means, if the response of the system due to load Px and Py is Rx and Ry respectively, then for the load (Px+Py), the response of the system will be (Rx+Ry).
(2) Time-Invariant means, the parameters of the system does not vary with time.

If the frequency response of the system is Time-Independent, then it suggests that the parameters of the systems are not changing over time, i.e., the system is time-invariant.

But it does not say anything about nonlinearity. Hence, in a nutshell, if the frequency response of a system is time-independent, this system is NOT an LTI?


How do I plot the spectrogram of a wav file? I am dabbing a little bit into signal processing on MATLAB. And I have the following question that I would like to know the answer to.

What has to be done : Plot the Spectrogram of an audio file.
What I know : Since its an audio file, I use wavread to read it and use the array returned to plot any of the graphs that I would like to. Only recently, I found that wav files store data on multiple channels. All the files that I have, have two channels. I know that the Spectrogram is supposed to give you an idea of the fundamental frequencies present in your audio file.
What I would like to know : What does the Spectrogram exactly do? What do the parameters mean? This is the structure of the Spectrogram function call : S=spectrogram(x,window,noverlap,nfft,fs). What do the window and noverlap parameters mean? How do I set these parameters? Are they input file dependent? After all this, how do I plot the Spectrogram of the audio file since it contains two channels? The function returns an error when I try to pass the two dimensional(since there are two channels) array returned by wavread. It says that I can pass only a row or a column at a time. How do I go about plotting the Spectrogram of my wav file?

Thanks a million, in advance! : )

This isn't a great question, but I'll try to answer:

  All the files that I have, have two channels.

A spectrogram can only show one channel at a time.  Either make one spectrogram for each channel, or average your channels together into one.  Here's an example of the syntax I found with a quick google search:

[y, fs]=wavread(fileName); % Read wave file
left=y(:,1); % Left channel
right=y(:,2); % Right channel

then you can do spectrogram(left)

  I know that the Spectrogram is supposed to give you an idea of the fundamental frequencies present in your audio file.

Not quite.  The term "fundamental frequency" means something different.

  What does the Spectrogram exactly do?

It shows the frequency spectrum of the signal and how it changes over time.  If the beginning of the signal is white noise (flat spectrum), and the end of the signal is a tone (spikey spectrum), the spectrogram will show how it changed from one spectrum to the other over time.  It does this by dividing the signal up into small chunks and calculating the spectrum of each chunk.

  What do the parameters mean? This is the structure of the Spectrogram function call : S=spectrogram(x,window,noverlap,nfft,fs). What do the window and noverlap parameters mean?

Did you read the documentation?

window is a window function that is applied to each chunk of the file before the frequency transform.  If you didn't use a window function, there would be abrupt changes at the beginning and end of each chunk, which produce wide frequency components that you don't want.  The window function is like a "fade-in, fade-out" envelope, so that you only see the spectrum of the middle of the chunk, and not the artifacts caused by chopping the chunk abruptly at the endpoints.

noverlap is the number of samples that each chunk overlaps the next.  The windowing function attenuates the signal to near 0 at the boundaries of each chunk, so if you didn't overlap the chunks at all, you would be emphasizing parts of the signal and ignoring other parts.  The overlap lets you measure the same part of the signal twice, once in the middle of the chunk where it's important, and once at the ends where it's mostly ignored.

  How do I set these parameters?

Basically you try different ones and see which is best at showing you what you want to see.  Different window functions are designed to do different things, so read about them and what they are meant to do.

  Are they input file dependent?

Only if you're looking for different types of things in different files.


FFT doesn't always have frequency in all frames I'm trying to do an FFT on live audio to see its spectrum.

So I have a sampling rate of 44100 Hz and from this take out 1024 bins on my FFT

So to test it I play a continuous tone at 15000 Hz and another one at 16000 Hz and look for it on my spectrum. I try to use higher ones to cut out noise from around. I have made the sound in a way there is no clipping

The thing is when I look at the individual frames of the FFT in each 2048 byte block of data I can see the peaks at these frequencies but some frames won't have the peaks, yet the others would.

How can this happen? How can I avoid it or what can I do to the way I'm reading the audio so that I can see the two tones in every frame?

You might be looking only at the cosine component (sometimes called the even or real component) of the FFT result, instead of the magnitude of each entire complex FFT result vector.  The real component will vary as your test signal goes in and out of phase with the cosine basis vectors (always even functions) of each frame, even if the complex bin magnitude stays constant (or varies only slightly due to windowing scalloping).

mag[i] = sqrt(real[i]*real[i] * imag[i]*imag[i]);


How can we calculate the frequency response of each filter of a Discrete Wavelet Decomposition filter bank? I am using the Multilevel Wavelet decomposition. The decomposition results in a filter bank such as the following:

Specifically, I am using one of the Daubechies wavelets, db7.

In reality Level 2 coefficients are the result of filtering by $\ g(n)$ then down-sampling and filtering by $\ h(n)$. Suppose that the whole operation I described above is a single filter ie. at each level the coefficients get generated by a single filter.
I would like to know how to calculate the frequency response of those filters at each level.

Is it correct to feed an impulse, i.e a Kronecker Delta function to the wavelet decomposition and then compute the Fourier Transform of the resulting coefficients ?

Thanks

PS: The same question has been asked here and the response was to get the convolution of the impulse response of $\ g(n)$ and an up-sampled version of $\ h(n)$. What puzzles me in that response though, is why divide by $\sqrt2$ the result of the convolution. Also, it seems like more code compared to just calling the wavelet decomposition routine with a Kronecker delta function.

The decision under the link you provided seems good enough. The main idea is:

we have impulse response of each filter
take the impulse response of the filter with the slowest sample rate as the input signal
route it back through the chain
perform spectrum estimation of the resulting signal

Division by $\sqrt2$ is for normalization only. For the given set of db coefficients their sum is $\sqrt2$, look here for explanation. It has no impact on the sense of the result I think.


gnuradio: How can I change the value of the tag exiting the "Correlate Access Code - Tag" block? Background: I have built a conventional demodulator in gnuradio and am trying to send payload data out of the flowgraph to an external Python script. I'm using gnuradio v3.7.12.

It seems to me that the optimal method would be to use the flowgraph to:

find the preamble and tag it
convert the data following the payload to a pdu
send the pdu as a message to the external Python code

I've built a flowgraph that successfully does this using the following blocks:
input bitstream->"Correlate Access Code - Tag Stream"->"Repack Bits"->"Tagged Stream to PDU"->"ZMQ PUSH Message Sink".

The problem is that the "Correlate Access Code - Tag Stream" block expects a very specific header between the preamble and the payload: 32 bits containing the payload size (16 bits repeated twice). This is fine if I am the one controlling the input signal. Any other signal I try to receive, however, is not going to have this header field.

So instead of the "Correlate Access Code - Tag Stream" block, I tried to use the "Correlate Access Code - Tag" block, but this block sets the length tag value to 0, which prevents the rest of the signal chain from working correctly. It seems like all I need to do is change the value of the length tag to my expected payload length, but unfortunately, after many hours of scouring the documentation regarding tags, messages and pdus, I have been unable to figure out how to do this.

Am I on the right track? I've attached my working flowgraph (with the "Tag Stream" block) as well as the non-working flowgraph (with the "Tag" block). I tried a number of things (Tagged Stream Align, Stream to Tagged Stream), but I couldn't get anything to work.

Working GRC: https://www.dropbox.com/s/ilv0ry1x4nvqelj/tag_stream_test.grc?dl=0

Non-working GRC (can't set length tag!):
https://www.dropbox.com/s/aoueiovuvj9ss8z/tag_test.grc?dl=0

Thanks for any help you can provide.
PC

I'm also hit by the same problem, and came up with "Tag Rewrite" python block to update the value.
https://gist.github.com/4a65be9c5139eb132f311b68808d36a6
However, even though I was able to update the length, "Repack Bits" block does not seem to pick it up. So while I was able to change the value of the tag, I still haven't accomplished my goal of taking the data out from the stream.


extended hamming code for the (7,4) Hamming code, I know how to decode a codeword using the syndrome decoding method, by just using the 3 parity bits, and it is correctable if there is just 1 error.

However, for the (8,4) extended Hamming code, I know how to decode the codewords correctly when there is no error or 1 error respectively. However, can we do that for the case of 2 errors?

Suppose the transmitted codeword is $\mathbf C$. If no errors occurred, or one error occurred, then $\mathbf C$ is the unique codeword
that is at distance $1$ or less from the received word $\mathbf R$, and the decoding algorithm produces $\mathbf C$ as its output. When there are two errors, there is not a
unique codeword that is at distance $2$ or less from the received word. One
of the codewords at distance $2$ is $\mathbf C$, but there is at least one other
such codeword. So,
the decoding algorithm must be devised to handle this case. The possible options
are

Declare a decoding failure, that is, tell the world that an uncorrectable
error pattern has been detected. A less shameful alternative is used in ARQ
(automatic-repeat-request) systems in which the receiver sends back a
NACK signal to the transmitter which then retransmits $\mathbf C$. Thus,
the outside world is not told that the transmission is being repeated
over and over till the decoder is able to successfully decode, though the
outside world may be able to infer this because it has to wait extra long
for the next data bits to come out of the decoder.
Choose one of the codewords nearest to $\mathbf R$ and produce this as
the decoder output. In this instance, the decoding might be correct because
the decoder fortuitously chose $\mathbf C$, or it might be incorrect.

It is important to remember that decoding algorithms can output a
wrong codeword even when they are fairly certain they got it right.
If $\mathbf C$ is transmitted and the channel error pattern $\mathbf E$
is such that $\mathbf C + \mathbf E = \mathbf R$ happens to be another
codeword $\hat{\mathbf C}$, then the decoder output will be
$\hat{\mathbf C}$ and the decoder might even declare that it is
confident that it has found the right codeword because it did
not find any errors in $\mathbf R$; it matched $\hat{\mathbf C}$
in every bit! Pride, as we are told, goeth before a fall.

Finally, turning to the OP's specific question about the
$(8,4)$ extended Hamming code, the decoder can
tell that an even number of errors have occurred if the
XOR sum of the $8$ received bits is $0$ and the rest of the
syndrome (the part that is exactly the same as in the
$(7,4)$ Hamming code)  is nonzero. So, whenever this
happens, the decoder can declare a decoding failure
or send a NACK to the transmitter, or output a codeword.
An easy way out in this last case
is to assume that the overall parity bit
is in error (and thus flip it) and then use the rest of
the nonzero syndrome to correct one error in the other
$7$ bits exactly as if they were from a $(7,4)$ (non-extended)
Hamming code. This will give a codeword $\hat{\mathbf C}$
that is at distance $2$ from $\mathbf R$ and it will even
be the transmitted codeword $\mathbf C$ if there were exactly
two transmission errors, one in the $7$ bits and one in the
overall parity
bit.


Filtering data so that only rising edge is left I have data that looks like this:

Sometimes the data has a higher point in the middle of the shallow slope

I want to find a way to filter the data such that it smooths it and leaves the first rising edge as the peak like so

Is a weighted average a good way to achieve this? or would that not do the trick? Are there better methods

Right now i have been using a weight average to smooth it but it leaves me with a peak in the center of the slope instead of towards the rising edge.

Try looking into posts on "onset detection", where typically "flux" (or first order difference, typically of spectral magnitude) is used to compute the "rate of change" of a signal. The high rate of change is the most salient feature of the images you drew above.

However, your third image seems to translate the maximum of the curve (to the left) is that your intent?


Using DFT Rotation for Frequency Down-Mixing Is it possible to perform Frequency Down-Mixing (sometimes referred to as Down-Conversion) by simply rotating the DFT sequence? If so what are the Advantages and Disadvantages of this method?

Is it possible to perform Frequency Down-Mixing ... by simply rotating the DFT sequence?

Yes.

If so what are the Advantages and Disadvantages of this method?

As compared to the standard method of just multiplying with a cosine or complex exponential, there are no advantages (hence, no one really does it). There are lot of disadvantages: it requires framing including overlap management and potentially windowing, it's computational expensive and calculates a lot of stuff that you actually don't need.


Haar transform definition I am trying to understand Haar transform of simple vector of numbers.

Some references I found say, that first level transform of $ X = (4,6,8,10,13,9,3,3)$ will be $\sqrt{2}(5,9,11,3, -1,-1,2,0)$. Others say, it will be $(5,9,11,3,-1,-1,2,0)$.

Which one is true? Or are they both correct in some way and am I missing something?

And important question, which one do I use when I am performing Haar transform on image and this vector represents one row/column of pixels in image?

The only difference between the two transform results is a gain factor.  The transform should have the same energy as the original data sequence, so we can use that to check for the right result.  The energy of $X$ is-

$$
\sum X^2 = 484
$$
$$
\sum (5,9,11,3,−1,−1,2,0)^2 = 242
$$
$$
\sum (\sqrt{2}(5,9,11,3,−1,−1,2,0))^2 = 2\sum (5,9,11,3,−1,−1,2,0)^2 = 484
$$
This shows that of the two possible transforms, $\sqrt{2}(5,9,11,3,−1,−1,2,0)$ has the right energy and is thus the right choice.


Implementing short-time Fourier transform I have a $512$ sample frame, lets call it X. The frame is further divided into four sub-frames of $128$ samples each. Each of these sub-frames have a different frequency.

For STFT, I use MATLAB's gausswin(N) function where N is equal to $128$ in my case to generate a Gaussian window, let's call it W. I zero-pad my windows W to make it equal to $512$ samples.

I now multiply my original frame X point by point with W to get y. I then take FFT of y to get Y. It works till here.

My question is, naturally the next logical step should be to shift my zero-padded window W to right by 128 samples and repeat the multiply and FFT procedure again.

But how do I accumulate these results?
Concretely, how do I convert it into a time-frequency plot?

admittedly, i have never used the buffer() function in matlab.  for doing STFT analysis and processing, i have just broken it into frames myself with code similar to this answer.  it's easier to see what's going on.

while you can use a gaussian window for analysis (i found it mathematically useful in this paper, email me if you want a pdf copy), you must use a complementary window, one that adds to 1 with the windows of its adjacent frames, for the reconstruction of the signal.  at least if you want this to pass the signal unchanged when the processing is set to "null".


A path to get the foundation of SP and DSP for a pure math I have a Ph.D. in pure math (interested in Harmonic analysis and operator theory). I am looking forward some proper references to lead me get the foundation of discrete/signal  processing more and more. Actually,  I had a  review of the  Heppenheim's books  (both signal and digital ones) and (rather) got what he is saying in these  (very nice) book. Now, I am going to develop my knowledge concerning this context and do need to cover  some more advanced ones.   Thanks in advances for you suggestions.
Probably based on my own field, I would  like to read on some texts whose approaches are focused  on  theoretical bases. However having some (proper) references which make me feel (realize) some real applications are in priority.

You have already read those Oppenheim's Signals & Systems, and Discrete-Time Signal Processing books.
I'm not sure what you mean by foundations but in some sense these two are also the foundations on signal processing. In other words, there are no (popular & successful) graduate level DSP books that discuss at an advanced level the same topics that are covered on them.
However the following books (or subjects) will enhance your understanding, or broaden your appreciation of the subject.
First of all, the very first graduate level course on DSP, Communications, and Control is called Linear System Theory which brings together all the undergraduate mathematical stuff from a new, advanced, deeper, and foundational point of view of the Hilbert (linear vector) Spaces, Linear Mappings, and Matrix theory. It does not have a definite book but a bunch of books on Linear Algebra & Matrices, Measure Theory, and Diferential Equations were used. Note that there's a control theory oriented bunch of Linear System Theory books (Desoer's crew) that I do not recommend for DSP, unless you will be designing control systems on the field. Signal processing does not make much use of state-space approach, unless it's absolutely necessary.
Then the second refresher/deepener is on Probability, Statistics and Random Processes. Fortunately it has two very strongly recomended books though:

Statistical Digital Signal Procesing _ Monson HAYES
Discrete RaNdom Signals and Statistical Signal Processing _ THERRIEN

The first book is a must read, the second one is following the style of Oppenheim series but is harder to follow, and less practical than the first.
Then the following books / subjects will be awating you :

Adaptive Filter Theory _ HAYKIN
Multiresolution Signal Decomposition _ AKANSU
Multirate Digital Signal Processing _ RABINER
Estimation & Detection Theory _ KAY
Pattern Classification _ DUDA
Theory and Applications of Digital Signal Processing _ RABINER & GOLD

Then the following applications will make your day :

Speech and Hearing for Communication - FLETCHER

Digital Processing of Speech _ Rabiner

Speech and Audio Signal Processing _ GOLD

Discrete-Time Processing of Speech _ PROAKIS

Advances in Speech Coding _ GERSHO

Two-Dimensional Signal and Image Processing _ LIM

Digital Image Processing _ GONZALES

Fundamentals of Image Processing _ JAIN

Signal Compression_JAYANT

Introduction to Data compression _ SAYOOD

Of course the list is by no means complete...


Building a spectral envelop of FFT'd audio I'm trying to figure out how to create a good spectral envelope of a signal.

Basically I am, at present, taking a windowed section of audio applying an FFT and generating a bar chart style representation of the magnitude of the spectrum.

However to add to this I'd really like to have a good spectral envelope overlaid over the top.  Unfortunately I just can't seem to come up with a good solution.

I tried a simple peak picking algorithm where I find every peak in the spectrum (I literally just look for points where the current value is greater than the values in the buckets either side).  This returns me a set of peaks and I then use a catmull-rom interpolation to draw the line over the peaks.

This follows the spectrum reasonably as you can see in the following screenshot:

However its not ideal.  In that image for example you can see how it doesn't follow the spectrum on the lefthand side.

When you zoom out:

You see that the envelope is pretty jaggy.  Its not terrible but is this really the best way to build a spectral envelope?

Is there a better way of picking the ideal peaks so that it truly envelopes the spectrum?  If so can anyone point me towards algorithms?

Also is catmull-rom interpolation the best way of doing it?  Can anyone suggest a better interpolation method and show me how to implement it?

The question is language agnostic but I'm writing in C++.

Edit: Ok I found an implementation of Auto-Regressive modelling.

I then implemented the method as follows:

    std::vector< float > coefficients( 3 );
    AutoRegression( &mSpecBuffer.front(), kFFTSizeDiv2, 3, &coefficients.front(), MAXENTROPY );

    mPeakBuffer.clear();

    mPeakBuffer.push_back( Peak( 0, mSpecBuffer[0] ) );
    mPeakBuffer.push_back( Peak( 1, mSpecBuffer[1] ) );
    mPeakBuffer.push_back( Peak( 2, mSpecBuffer[2] ) );

    int x       = 3;
    int xMax    = kFFTSizeDiv2;
    while( x < xMax )
    {
        const float k3  = (mPeakBuffer.end() - 3)->peakHeight;
        const float k2  = (mPeakBuffer.end() - 2)->peakHeight;
        const float k1  = (mPeakBuffer.end() - 1)->peakHeight;

        const float k   = (k1 * coefficients[0]) + (k2 * coefficients[1]) + (k3 * coefficients[2]);

        mPeakBuffer.push_back( Peak( x, k ) );
        x++;
    }

This gives me the following result:

Which to my eye looks a HELL of a lot better.

Zoomed right out it still looks pretty damned good:

So I just want to check, if this is correct?  If so I'll work along this path further :)

Two techniques commonly used for recovering the spectral envelope of an audio signal:

Low-order AR modeling. Estimate your AR coefficients, then plot the magnitude of the frequency response of the estimated all-pole filter.
Cepstral analysis. Take the logarithm of the magnitude of the FFT output. Take the FFT, truncate to keep only the first coefficients, take the inverse FFT, take the exponential. This is equivalent to applying a low-pass filter on the spectrum (with a log-scale magnitude).


An explanation of parts of the MUSIC Algorithm for someone without a background in signal processing? I've been looking at every tutorial and paper I can find, but I keep losing the intuition behind what's going on when they get into the matrix operations.

There is a step in which the signal vector is multiplied by its hermitian transpose and the correlation matrix of the product is computed. What do both of these steps accomplish from a high level?

Where is the 'e' term arising from in the steering vector?

I also don't understand the last part about how the orthogonality of the signal and noise subspaces is used to disentangle signals. We use a function in which we have 1 divided by the product of vectors representing the signal and noise subspaces. If they equal 0, wouldn't this function be undefined since it would be equal to 1/0?

First of, if you can read German: I was author on a thing, back in the day.

There is a step in which the signal vector is multiplied by its hermitian transpose and the correlation matrix of the product is computed. What do both of these steps accomplish from a high level?

So, first of all, MUSIC can be used for a couple of things. Originally, it was used for direction finding (underwater microphones, if I remember correctly), and later was used for spectrum estimation. The idea is this:
You have two inputs. In the direction finding usage, those are two different microphones, in the frequency estimation case, it is but one signal.
You then estimate the covariance matrix – I'm going to focus on the frequency estimation part for now, so you're focussing on the autocovariance matrix – which basically says "If I shift my signal vector by N samples, this is how similar it looks" ("similarity" = dot product in each element of the matrix). That's your dyadic produce (ie. multiplication with conjugate transpose) – nothing but an estimator for the autocovariance matrix.
Then you go

"OK, let's abstract this and take it from a pure linear algebra perspective: these column vectors of this $M\times M$ matrix describe my signal, and knowing that they must be linearly independent (noise being uncorrelated, and present), they form a base of the receive signal space.
What about we decompose that space into one subspace with only noise vectors, and one containing only signal vectors?
If we do that, the base vectors of the signal subspace must be the periodicities contained in my signal.

So you go ahead and do that with an appropriate Eigenvalue decomposition, sort the eigenvectors by magnitude, draw an arbitrary line between noise and signal, and say "ok, the eigenvectors belonging to eigenvectors above that threshold are signal, the others noise representatives".


Modification of Accelerometer Data for Feature Extraction I have accelerometer [x,y,z] data recorded from various activities.  I am not looking at any filtering techniques (which I have seen within this forum) because I want to be able to measure the amplitude of the noise for each activity sample set using a global time-domain feature average loudness.

I am curious to know if I can center the signal around 0g, similarly to the X Data Minus Mean vector plot, but without the x-axis orientation i.e. the signal shifting up and down?

I have the x-axis data from one of my activities below:

I was hoping that by subtracting the mean from the x-axis data would yield this, but all it has done is shifted the entire signal down to around 0g.

The reason I ask is because I am unsure if computing a time-domain feature such as average loudness will work with the signal shifting up and down?

Average Loudness/Energy:

Your data is composed of a majoriy of low-activity regions, with higher-amplitiude  sparse/concentrated activities. The bad news is that most of the methods to achieve your goal are likely to somehow filter your data. The good news is that non-linear filtering could help you.

If you are not concerned with calculations and real-time yet, I would suggest first a centered median filter to extract a robust $0 g$ estimate, and remove it from the data. In other words: choose an appropriate left-right span of $K$. In each window around index $k$, computes the signal $\hat{s}(k) = s(k) - \textrm{median} [s_{k-K},\ldots,s_{k+K}]$.

A median is a non-linear filter minimizing the $ell_1$ norm (like the $ell_2$ norm minimization yields the standard mean). It is somewhat robust to outliers, and since yours are both above and below your thought $0 g$ estimate, a running median (computed around each sample) could extract this 0-level reference.

Looking at your plots, start with $K=50,100,200$ to see if this starts to provides you with what you are looking for. Then, more involded techniques could be considered.


Find Difference Equation From Frequency Response Given the frequency response $H (jω) = csc(ω)$ how would I go about finding the difference equation?  I know I can find $h[n]$ from here but after that I feel like I'm missing something.

so, if we're looking around for a discrete-time difference equation, let's modify the notation slightly to make it look like it's a frequency response for a discrete-time system:

$$  \begin{align}
H \left( e^{j \omega} \right) & = \csc(\omega) \\
& = \frac{1}{\sin(\omega)} \\
& = \frac{2j}{e^{j \omega} - e^{-j \omega}} \\
& = \frac{2j \ e^{-j \omega}}{1 - (e^{-j \omega})^2} \\
\end{align}$$

so i guess that means that
$$H(z) = \frac{2j \ z^{-1}}{1 - z^{-2}}$$

so can you get a difference equation from that?


Non-coherent discriminator I'm reading a paper about designing a non coherent discriminator for GNSS receiver. I came to the equation

I'm not sure about that symbol <
how would I write this in a small pseudocode ? for example C++/python ?

from the paper

When considering PLL with extended integration time, several correlator outputs, obtained from subsequent portions of the input signal, are evaluated and used for producing an improved phase estimate. For this reason, the index k = 0, 1, ..., K − 1 is introduced and the different correlator outputs are denoted as Pk = PI,k + jPQ,k = dkAc exp{ jφ} + ηP,k.

also I'm wondering how would I describe the steps for that loop archietecure in pseudo code ?

the loop architecture depicted in Fig. 2 is obtained. The input signal is, at first, integrated over one bit period; the complex correlator outputs, Pk, are then squared, further integrated, and fed into a four-quadrant arctangent that provides the discriminator output.

That $\angle$ is "angle". As the text right below it indicates, it's the phase of a complex number, in this case the phase of the sum of squared $p_k$.

Your second question: I really don't think pseudocode makes this any clearer, but sure:

factor = 1
sequentially for each timestep n:
    s = r[n] * factor
    t = integrate_over_last_N_and_decimate_by_N(s)
    u = remove_bits(t)
    v = non_coherently_integrate_over_last_K(u)
    w = non_linear_discriminator(v)
    x = loop_filter(w)
    y = nco(x)
    factor = generate_carrier(y)

I mean, that's a pretty straightforward block diagram.


How to keep the same power (sonority) between pure sines? I have a list of sine waves in 1/3 Octave (63,80,100...1k,1.25k,etc) that I want to reproduce with the same power. Right now I just have a global gain that applies to all this freqs, but of course the high frequencies are louder that the lower ones.

I am not sure how to keep the same power between them, I suppose I can use the principle of the pink noise and convert the gain (normalized 0-1) to dB, lower 3dB per octave as I go up in the list and go gack to normalized. But I am not sure if this is the correct way to achieve what I need.

How can I Keep the same power/Sonority?

Edit:
By equal power I mean they sound with "the same volume" to person listen to them.

For info, the user will be using mostly headphones. This is running on mobile phones.

From the comments on the question I gather that you want the individual sine waves to appear to be at the same volume. Unfortunately, there is no single EQ setting you could apply to achieve this. As for the reasons, let's look at the equal-loudness contours:

Source

What these contours tell us is that for a specific level of phones there exists an EQ setting that would make your sine waves sound at equal loudness. Unfortunately, the human perception of loudness not only depends on frequency but also on SPL—what that means is that your EQ settings that you meticulously set up for a certain level of phones won't work any more if you change the overall gain of your playback system. To compensate for the SPL-dependency of the human hearing you would need an adaptive EQ that takes into consideration the SPL of the signal at 1000 Hz (or any other reference frequency) by measureing it with a sound level meter and setting the EQ parameters accordingly.


Digital AGC signal level Analog AGC is mainly used to increase the signal level to match the ADC sensitivity? Is Analog AGC also used to attenuate the signal?.
Digital AGC is the first block used after ADC. Is the main purpose of digital AGC is to provide constant level to the signal processing receiver chain( RRC, CIC,  demod....)
So analog AGC output is range of signal levels but digital AGC always gives out a fixed level?

Analog AGC is mainly used to increase the signal level to match the
ADC sensitivity? Is Analog AGC also used to attenuate the signal?.

For "Analog AGC" we use analog control components (typically voltage variable amplifiers and voltage variable attenuators) to adjust the receiver gain, and we use an analog power detector to measure the total band power (within the bandwidth of the filtering up to the detector). This measured power is compared to a target which generates an error signal (often with an op-amp in an analog circuit) and the error is integrated in a loop filter (also often with the same op-amp) which then integrates to the control voltage needed to make the error go to zero (and the loop filter can be further modified for higher order loop control). In the end the band power AS MEASURED AT THE DETECTOR will be constant except for dynamic conditions during power transitions.  Depending on what is used in the receiver for gain control, the AGC can amplify or attenuate all the signals in band as allowed.

Digital AGC is the first block used after ADC. Is the main purpose of
digital AGC is to provide constant level to the signal processing
receiver chain( RRC, CIC, demod....)

Digital AGC is not necessarily the first block after the ADC, or at least the level detector part of the AGC loop. A significant benefit of doing a subsequent AGC digitally is that you can measure the level (which sets the level of the overall waveform within the bandwidth at the detector) after more thorough filtering of the channel signal of interest and rejecting out of band interference (which could include other channels in a multiband receiver). This interference would reduce the signal of interest in band when AGC'd (power measured) earlier in the chain. Consider the above case earlier for the analog AGC where we have a very strong interference signal near the channel bandwidth of interest, so close that filtering in the analog would be costly. Assume too that we have a high dynamic range ADC with excess dynamic range allowing us to shift some of this filtering challenge to the digital. The analog AGC's purpose is to avoid excessive clipping at the ADC input (or other similar level sensitive components in the receiver chain), so with the strong interference our signal of interest will effectively be significantly attenuated, give the interference will set the full scale level at the ADC input.
So the Digital AGC's purpose is to further bring this signal back to the proper full scale level (and minimize subsequent receiver processing) by selecting the bandwidth of interest and removing the remaining interference digitally, and then once removed, bring the channel signal back to the proper full scale level (some determined back-off to minimize clipping in the waveform but also reduce the subsequent operations required assuming a fixed point design, meaning minimize datapath requirements).
Also we must be careful to consider all components between the adjusted signal level and the detector (assuming filtering prior to the detector means that the actual power in the signal affecting the components in line can be much higher leading to saturation, clipping and gain compression.

So analog AGC output is range of signal levels but digital AGC always
gives out a fixed level?

Not quite. Both analog and digital AGC provide a fixed level of channel power, as measured by the relevant power/signal detector used. What may more likely occur in the analog is strong interference controlling the AGC, so the AGC sets its level based on the interference and not our signal of interest. The AGC only measures power and does not recognize signal from interference: we help that with filter design.
There are also systems where the level is measured digitally and the front-end gain is controlled but it is important to monitor the total signal level at the input the ADC to optimize the trade between quantization noise and clipping, maximizing the SNR of the signal with the proper AGC set-point at the ADC input. (Some clipping is often actually better than backing off too far and introducing more quantization noise in the process).
The AGC loop consists of a detector and a way to adjust the power along with an appropriate loop filter as a control loop. The power adjustment can be done as gain or attenuation or both, depending on what is available in the design. The detection and gain adjust and loop filter can all be done digitally or analog or both and as suggested above it can be a multi-loop design, and the gain can also include stepped gain control.
For more info see these posts: How to adjust receiver gains to avoid saturation and quantization noise to optimise post digital processing?
Role of Power Metric and LPF in AGC(Automatic Gain Controller)
AGC noise level


"cascaded" cross-correlation so, first let's do this for the circular correlation case (it might be easier).

let's say we have three real and periodic functions or signals $x(t)$, $y(t)$, $z(t)$ all with period $2 \pi$:

$$ x(t+2\pi) = x(t) \quad \forall t $$
$$ y(t+2\pi) = y(t) \quad \forall t $$
$$ z(t+2\pi) = z(t) \quad \forall t $$

define the cross-correlation between any two, say $x(t)$ and $y(t)$ as

$$ R_{xy}(\tau) \triangleq \frac{1}{2\pi} \int\limits_{-\pi}^{\pi} x(t) y(t+\tau) \ dt $$

let's say that the maximum value of $R_{xy}(\tau)$ occurs at $\tau=\tau_{xy}$

$$ R_{xy}(\tau_{xy}) > R_{xy}(\tau) \quad \forall \tau \in (\tau_{xy},\tau_{xy}+2\pi) $$

(since $x(t)$, $y(t)$, and $z(t)$ are all periodic with period $2 \pi$, so also are any cross-correlation between any pair.)

$$ R_{xy}(\tau_{xy}) = R_{xy}(\tau_{xy} + k \ 2\pi) \quad \forall k \in \mathbb{Z} $$

so here's the question: except for a difference of some multiple of $2\pi$, is

$$ \tau_{xz} = \tau_{xy} + \tau_{yz}$$

or

$$ \tau_{xz} = \tau_{xy} + \tau_{yz} + k \ 2\pi $$ for some integer $k$?

i think this should be true, but it would be nice if it is proved to be so.

I'm afraid your statement isn't true. This can best be seen in a suitable choice of basis, one that simplifies the cross correlation. This basis if of course the shift invariant periodic Fourier basis on your support interval.

Let's label the basis vectors $F_n$ for integer $n$. The cross correlation of two different such basis vectors vanishes, because the vectors are orthogonal and shift invariant: $\langle F_n,F_m \rangle =0$ for $m\neq n$. If we slightly generalize the cross correlation to work for complex signals as in $$\langle A,B \rangle (\tau):=\int_{-\pi}^\pi A^*(t)B(t+\tau)dt$$ and assume normalization of the Fourier basis, then $\langle F_n,F_n \rangle (\tau)=F_n(\tau)$.

With this, we can construct your $x(t)$,$y(t)$ and $z(t)$ signals like follows:
$$x(t) := c_{xy} F_{n}(t) + c_{xz} F_{k}(t)$$
$$y(t) := c_{yx} F_{n}(t) + c_{yz} F_{m}(t)$$
$$z(t) := c_{zy} F_{m}(t) + c_{zx} F_{k}(t)$$

with mutually different integer indices $k,m,n$. The complex coefficients $c_{ab}$ for $a,b\in\{x,y,z\}$ can be chosen arbitrarily and independently.

The corresponding correlation functions are:

$$\langle x,y \rangle (\tau) = c^*_{xy} c_{yx} F_n(\tau)$$
$$\langle y,z \rangle (\tau) = c^*_{yz} c_{zy} F_m(\tau)$$
$$\langle x,z \rangle (\tau) = c^*_{xz} c_{zx} F_k(\tau)$$

So each cross correlation function is a complex sinusoid with a phase factor in front. The phase factors and the sinusoids of all three correlation functions can be chosen absolutely independently. If you want to talk about maxima instead of phase factors you can take the whole setup to the reals by either adding the conjugate basis functions or just take the real part of these correlation functions. Either way, you can chose the coefficients to give you any set of maximum lags $\tau_{xy,zy,xz}$ without the requirement of the connection between them you conjectured.

So as requested let's get more explicit. The construction principle that guarantees independence of the three cross correlation functions should be relatively obvious: The orthogonality of the Fourier basis allows us to construct the signals so that the signals $x,y,z$ only overlap pairwise. In other words, the example above demonstrates that there is no transitivity for cross correlations that would allow general statements about $\langle x,z \rangle$ from $\langle x,y \rangle$ and $\langle y,z \rangle$.

The simple construction from above makes them overlap with just one basis function, which is a complex sinusoid. This is probably why the result cannot be immediately seen to apply to the original question. If you want maxima we need to use real numbers, and single real basis functions also do not have a single maximum, so we need to have more than one basis function of overlap between the signals (except for $F_1$, which already only has one maximum per period).

The simplest full example that I can come up with right now and that follows the design principle described above is:
$$x(t)=\sin(t)+\cos(4t)+\cos(5t)$$
$$y(t)=\sin(t)+\sin(2t)+\sin(3t)$$
$$z(t)=\sin(2t)+\sin(3t)+\sin(4t)+\sin(5t)$$

The base indices that overlap between the signals are $1$ for $xy$, $2$ and $3$ for $yz$, and finally $4$ and $5$ for $xz$. The latter have a phase shift introduced to realise your requirement. The overlaps group adjacent basis functions so that we have a single guaranteed maximum.

At this point it's probably easier to inspect the resulting correlation functions graphically.

First, $\langle x,y \rangle$:

$\langle y,z \rangle$:

and finally the independent third pair $\langle x,z \rangle$:

The maxima of the first two cross correlations are clearly at $\tau=0$, while the third pair correlates best at some $\tau>0$.

This example works in both discrete and continuous time, as it only uses properties of the Fourier basis that are realised in both domains.


Average power of finite length discrete signal I know that the power of a discrete signal is given by:

$$P(x) = \lim_{N \rightarrow \infty}\frac{1}{2N + 1}\sum ^N _{n = -N} | x[n]|^2$$

However, I would like to know whether the power of a discrite time finite length signal is defined as:

$$P(x) = \frac{1}{N - 1}\sum ^N _{n = 1} | x[n]|^2$$

Or as,

$$P(x) = \frac{1}{N}\sum ^{N - 1} _{n = 0} | x[n]|^2$$

Being $N$ the length of the signal, the number of samples of it.

A finite length signal has a power of zero. This becomes obvious if you just use the first formula in your question and apply it to a signal of finite length. Finite length signals are so-called energy signals, and their energy is given by

$$E_x=\sum_{n=0}^{N-1}\big|x[n]\big|^2\tag{1}$$

where $N$ is the length of the signal (and we assume that the signal is defined in the index range $n\in[0,N-1]$).

The expression

$$P_x=\frac{1}{N}\sum_{n=0}^{N-1}\big|x[n]\big|^2\tag{2}$$

is the power of an $N$-periodic signal, as explained in this answer.


Create a color image in Matlab with sinusoidal patterns I want to create a color image in Matlab whose components have sinusoidal patterns (okay even if they have different frequencies).

To create a black and white image I did:

fs = 0.08;
W = 256;
H = 128;

img1 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img1(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

colormap(gray(256));
imshow(img1);
title('img1');

Simply do it for each color channel, I mean for Red,Green and blue, for each form a separate matrix and then concatenate those, to form a 3D color image.

fs = 0.08;   //FOR RED CHANNEL
W = 256;
H = 128;

img1 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img1(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

fs = 0.16;    //FOR Green CHANNEL

img2 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img2(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

fs = 0.01;    //FOR Blue CHANNEL

img3 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img3(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

Color=zeros(128,256,3);
Color(:,:,1)=img1;
Color(:,:,2)=img2;
Color(:,:,3)=img3;
imshow(Color,[]);


Are there theories to manage states of signal volatility around a limit Assume input signal and discrete output. The signal can cross the limit at 1000 times per second, but it is desired that the discrete output changes at the most once every 15 seconds.  See image below:

$U$ is the Signal, $A$ is the comparator output, $B$ is the filtered output.  The goal is to filter out the output changes in region $R$.

Are there theories to dealing with volatility around the limit such as in region $R$ ?
At this point I am thinking of an algorithm where state is only changed if the signal is persistent for time period $T$, is there a name for such an algorithm?

In electronics, a Schmitt trigger converts an analog signal to a digital signal by using two different thresholds (limits as you call them), one for an output state change from 0 to 1 and another for a change from 1 to 0. The first is somewhat higher so that if the input fluctuates between those two thresholds, then there will be no change in the output. The general concept of reluctancy to change state is called hysteresis.

It can also help to lowpass filter the signal first, which is often used in debouncing a mechanical switch, to remove fluctuation very much like in your picture.


Inverse Z-transform with complex conjugate poles I was computing an inverse z-transform here, and I am facing some problems.

So, the z-transform is:

$$ X(z) = \frac{2+3z^{-1}}{1 - z^{-1} + 0.81z^{-2}} , |z| > 0.9 $$

I found the following poles:

$$ p_1 = 0,5 + j0,7483 $$
$$ p_2 = 0,5 - j0,7483 $$

and then, the following residues:

$$ R_1 = 1 - j2,6727 $$
$$ R_2 = 1 + j2,6727 $$

Rewriting $ X(z) $ with the poles in polar coordinates:

$$ X(z) = \frac{1 - j2,6727}{1 - 0,9 \angle 56,25^{\circ}} + \frac{1 + j2,6727}{1 - 0,9 \angle -56,25^{\circ}}  $$

Through the magnitude of the poles, I can see that they are both inside the ROC, then, I can write $ x(n) $ as:

$$ x(n) = (1 - j2,6727)(0,9 \angle 56,25^{\circ})^{n}u(n) + (1 + j2,6727)(0,9 \angle -56,25^{\circ})^{n}u(n) $$

I know that it is possible to write $ x(n) $ in a form without complex numbers. I know this, because it was asked and I was able to plot the figure in Matlab using the following code:

delta = zeros(30, 1);
delta(1) = 1;

a = [1 -1 0.81];
b = [2 3];

x = filter(b, a, delta);
n = 0:29;

figure
stem(n, x);

My question is, what I did up to now is right? If it is, How can I compute $ x(n) $ without complex numbers? Is there any tips that you can give me to follow?

Thanks in advance.

Since the given $\mathcal{Z}$-transform

$$X(z)=\frac{2+3z^{-1}}{1-z^{-1}+0.81z^{-2}},\quad |z|>0.9\tag{1}$$

has real-valued coefficients, its poles and zeros must be either real-valued or they must occur in complex conjugate pairs. The poles of $(1)$ are a complex conjugate pair:

$$p_1=0.5 + 0.74833j=re^{j\phi}=p\\
p_2=0.5 - 0.74833j=re^{-j\phi}=p^*$$

with $r=0.9$ and $\phi=0.98177$.

You can proceed from the partial fraction expansion of $X(z)$, which you correctly obtained:

$$X(z)=\frac{A}{1-pz^{-1}}+\frac{A^*}{1-p^*z^{-1}}$$

With the given region of convergence (ROC) $|z|>0.9$, we know that we're looking for a causal sequence $x[n]$:

$$x[n]=\left(Ap^n+A^*(p^*)^n\right) u[n]\tag{3}$$

With $p=re^{j\phi}$ and $A=|A|e^{j\phi_A}$, $(3)$ can be rewritten as

$$x[n]=|A|r^n\left(e^{j(n\phi+\phi_A)}+e^{-j(n\phi+\phi_A)}\right)u[n]=
2|A|r^n\cos(n\phi+\phi_A)u[n]\tag{4}$$

which is of course real-valued.

With the given values of $p$ and $A$ the impulse response looks like this:


How to transform large set of audio samples using FFT? Well, I am new to DSP so I apologize for the basic question. I am trying to transform audio 16bit PCM (which I store in normalized 0.0-1.0 range) samples into frequency domain using fft.

I understand that FFT transforms N points of time domain into two N sized frequency domain array. My problem is determining the size of N. For example if I have a 5mins PCM signal sampled at 8000 Hz/sec I can't find a way to get the spectrum of the whole signal, and not several arrays of size (e.g. 1024, 2048 etc..).

So my question is:

Can I transform the whole signal at once ?
If I cannot, is there away to combine the spectrum of different N sized segments?
Is it reasonable to study the spectrum of a whole large signal at once ?

Can I transform the whole signal at once ?

Yes, this would be a 2400000 points FFT - something that is easy for a modern desktop computer (note that the sample rate is measured in Hz or samples/sec, not in Hz/s).

  If I cannot, is there away to combine the spectrum of different N sized segments?

Averaging the spectrum obtained on successive short windows is a common spectral estimation technique - but I'm not sure this is really what you need here.

  Is it reasonable to study the spectrum of a whole large signal at once ?

It makes sense if the signal is stationary. I can think of some applications in which this would be the case (for example you record 5 min of continuous hum/noise generated by a machine), but this would be an exception more than a rule - and I can't think of a practical reason why one would have to perform a FFT on such a large chunk of signal.

For signals whose characteristics change with time - such as speech or music - spectral analysis is more commonly performed on short segments. It all depends on what you want to do with the result of the FFT. If it is to extract a property of the signal (such as its pitch, or a metric like spectral centroid), make sure it is reasonable to assume that the property you want to extract is constant over the duration of the signal. If not, break down the signal into short segments during which this property can be considered constant.


Creating an Eye Diagram for QPSK Simulation I am trying to implement an eye diagram for an application, where the input signal is QPSK. However, I feel that there is some fundamental concept
concerning these plots that I am missing. Several definitions and descriptions for these diagrams that I have seen are all along the same line:

  The eye diagram repeatedly overlays the time width of n symbols

Which sounds straight forward enough (though there are possibly variations?), but I'm not sure that is all there is to it.

The simulated input signal that I am testing with:

modulated QPSK, generated from random symbols
No raised cosine, RRC, or any filtering
No noise added

Eventually, I will modify the signal (such as adding filtering) to see the effects on the system.

This is the image I get when overlaying the symbols in time (showing 1 symbol per trace):

Since I've never seen any example that looks like this, I tried looking at other variations. Using the eyediagram function in Octave, it produces the barn door (breaking the signal into
the real and complex):

The points on the "door" are just the constellation points received (not the samples). So when I see examples like this, with a noisy signal:

I don't have enough reputation to post another image here, though it would help to explain my question. Similar to the above image: for a noisy signal, the lines of the diagram are fuzzy. It would indicate that there are more analog points used to fill in the plot

Where are the other points coming from to create the noise, if not the received waveform? Then there are the horizontal components of the image. How is it possible to get the horizontal components from a QPSK waveform? Even if separated into I and Q representations? Again, it makes sense when connecting the received constellation point at T intervals, but I do not see how to get this when plotting the signal itself.

What am I missing or not understanding here?

EDIT

I updated the diagram to plot only the received symbols. Previously I was plotting the received
waveform, which was a modulated signal (that is how the first image was produced). Below are two
diagrams showing just the in-phase plot. The second one has noise added:

Without noise, the lines are straight as PSK mentioned in the comments. With noise, the lines are
still straight. Which is the other part of the question. Looking at an example here:

(this is from Matlab example at https://www.mathworks.com/help/comm/gs/scatter-plot-and-eye-diagram-with-matlab-functions.html)

Where do the smooth transitions come from? The lines are not straight. There are other in-between points that are filling in the diagram. Where are they coming from?

Your eye diagrams for the inphase and quadrature baseband signals are perfectly correct. For QPSK with rectangular signals and perfect matched filtering, the inphase and quadrature signals are indeed just BPSK with rectangular signals, and the matched filter outputs can be deduced from, for example, the figure (reproduced below) at the end of this answer of mine.

If you sketch the matched filter output between $T$ and $5T$, overlaying the signal in chunks of length $2T$ centered at $2T$ and $4T$, you will see that the signal is flat at the lower level when two $1$s are transmitted in succession. Similarly, the signal will be flat at the upper level if two $0$'s are transmitted in succession. For transitions between two data bits, the matched filter output crosses over.

As to where the smooth transitions instead of straight ramps come from, if the modulating signal is not a rectangular pulse, or if the matched filter is not quite matched perfectly, the filter output can well be smoothly transitioning from one level to the other. For example, if the BPSK pulse were one cycle of a sinusoid of period $T$, specifically
$$s(t) = \begin{cases}\sin\left(\frac{2\pi t}{T}\right), &0 \leq t < T,\\
0, &\text{otherwise,}\end{cases}$$
then the matched filter output would be a perfect sinusoid of period $T$ between two successive sampling instants if the two data bits were the same, and make a smooth transition from peaking at $nT$ to peaking with opposite polarity at time $(n+1)T$.
See, for example, the figure below.


Half Subcarrier spacing in LTE Uplink First post here, trying to make it as good as possible. Do not know how to make the Python code look nice - any tips? :) Anywho - to the main problem:

I am currently working on simulating the uplink of LTE starting with the SC-FDM(A) modulator. Now, the "first edition" is the OFDM modulator implemented in Python and basically:

import numpy as np

import numpy.fft as npf

numSymbols = 12*10 # 10 resource blocks, each containing 12 subcarriers.   FFTsize = 1024

tmp1 = np.array(np.random.randint(0,2,numSymbols))*2-1   tmp2 = np.array(np.random.randint(0,2,numSymbols))*2-1

inputSymbols = (np.array(tmp1) + 1j*np.array(tmp2))/np.sqrt(2) # QPSK modulation

TxSamples = npf.ifft(inputSymbols,FFTsize) # Zero-padded IFFT.

TxSamples is the OFDM modulated signal of some random binary input. The subcarrier mapping is implicit in the IFFT due to the zero-padding.

My issue is the fact that LTE Uplink uses half a subcarrier shift in relation to DC. My implementation is build upon the assumption that each entry in the IFFT is spaced with 15kHz (the required subcarrier spacing according to LTE) and this is also the case for two other implementations I have found here and here.

So, how can this half subcarrier shift (7.5kHz) be incoorperated into the simulation?

Let me know if you need more information - do not really know what is to much or to little information regarding my issue.

Thanks in advance.

** UPDATE: **

So, I tried to do the frequency oversampling by simply adding a zero in between samples. This resulted in me having a replica in the time domain. Cutting away this replica got me back to the starting point.

What I have done in order to solve the problem is basically mathematical driven. From the formula for the DFT it is possible to derive that you can simply multiply the input symbols by a complex exponential. This will shift the frequency. The Python code is basically (overwriting the last line in the example above):

ism = []
for m in range(len(inputSymbols)):
    ism += [inputSymbols[m]*np.exp(1j*np.pi*m/FFTsize)]
TxSamples = npf.ifft(ism,FFTsize)

which yields the correct time domain representation, with no sorta explicit oversampling. To verify scipy can be used as:

import scipy.signal as sps
import pylab as pl

ts1 = sps.resample(TxSamples,len(TxSamples)*2)
pl.figure(1)
pl.plot(10*np.log10(abs(npf.fft(ts1))))

which should show that the signal has now been oversampled by two and that the frequency bins are located correctly, while the time domain representation is not broken.

I suppose this way one makes sure that both domains are "the same". Since by simply oversampling in the frequency domain it is not possible to represent the time domain signal correctly. I do not know - it seems that this "hacked" solution works. However kind corrections are always welcome.

For implementing the frequency shift of 7.5 kHz you will have to apply oversampling at some point in your system. There are two principal methods to achieve this: in frequency domain or in time domain.

Frequency domain

You could simply apply an additional zero padding in frequency domain by using a 2048 IFFT and by setting every second subcarrier to zero. You then shift all subcarriers by one frequency bin to achieve the frequency shift. This will maintain the frequency spacing of 15 kHz but the bandwidth of each subcarrier will be divided by two compared to the 1024 IFFT.

Time domain

You apply two-times oversampling to the time domain signal at the output of the IFFT. This can be done by padding a zero sample after each sample followed by a lowpass filter or some other interpolation method. The frequency shift can then be achieved by multiplying the resulting signal by $\exp(j\Omega_1 n)$, where $\Omega_1$ is the 7.5 kHz normalized to the new sampling frequency (double the original sampling frequency) and $n$ is the discrete time.

Which method to choose depends on how close your simulation should be to the hardware implemenation. If the system you model has some intermediate frequency oscillator to do the frequency shift, the time domain method can model this. But if your digital-to-analog converter is fast enough and if you can implement a 2048-IFFT at the required speed, this is probably the way to go and you should therefore model it accordingly.


What is the meaning of the DTFT of the unit impulse sequence? In an exercice, I'm asked to draw the $X_{imp}(\omega)$ Discrete-Time Fourier Transform (DTFT) of the $x_{imp}(n)$ unit impulse sequence defined as:
$$
x_{imp}(n) = \begin{cases}
  1 & \text{if } n = 0 \\
  0 & \text{otherwise}
\end{cases}
$$

Thanks to another question, I think I understand the difference between the DFT and the DTFT.

So, given $X(\omega) = \sum_{-\infty}^{\infty} x(n)\cdot{\rm e}^{-j\omega n}$, I end up with $X_{imp}(\omega) = 1$ because all terms of the summation would be $0$, except at $n=0$ where $x(n)\cdot{\rm e}^{-j\omega n} = 1$.

If my answer is correct, does that mean the unit impulse sequence would have the same amount of "energy" (← is this the right word?) for all frequencies, including audio, radio frequencies X-Rays, light or whatever else. That seems somewhat irrealistic to me. So, where is the flaw in my reasoning?

Your calculation is correct.

Your signal has its energy equally distributed in all frequencies.

But note that in discrete-time, "all frequencies" means from $\omega = -\pi$ to $\omega = +\pi$. Frequencies beyond that range exist but result in the same signal as one with a frequency within the range.

For example, $\cos(0\times n) = \cos(2\pi \times n)$. Similarly $\cos(0.1 \pi \times n) = \cos(2.1 \pi \times n)$, and so on.

So the energy of the signal is calculated integrating in a range of $2\pi$.

The formal result is known as Parseval's theorem, which states that energy measured in time domain is equal to energy measured in frequency domain:

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

Note: this is a particular version of Parseval. The more general version states that the projection (think cross product) of one vector onto another is the same calculated in time domain or frequency domain, which is reasonable since Fourier can be seen as a change of base, and a change of base should not alter projections.


2d dft of a scattered signal I want to take the DFT of the data shown in the image. Note: my data consists of $(x,y,z)$ values for which I have shown $(x,y)$ in the image (i.e. their location). '$z$' is the amplitude of the data points which I haven't shown. I have some questions:

Is 2d DFT for this data the same as taking multiple 1d DFT for all the rows? (Should I venture into the world of 2d Fourier to analyze this data of mine or is 1d analysis sufficient?)
As it can be seen ,my data is highly scattered, do I need to interpolate it before applying DFT? If yes, which interpolation would you suggest ? Or, can I use a non uniform DFT ?

I am a newbie to DSP and I understand my questions are elementary. So if you feel it does not need to be answered here, you may direct me to a web-page or a link to a question-answer/discussion which addresses this issue.

Also no data does not exist between the blue dots i.e the blue dots represents the points where I have the height (amplitude) value and the in between 'white' points I do not have the height (amplitude).

Without knowing what you want to get, I will assume you are interested in periodicities of the elevated dot patterns. They can give you access to roughness measures, texture attributes.

Visually, 2D periodicities are apparent from your data point $(x,y)$ locations. At different scales though: beside a potential periodic base pattern, little streaks of 5 to 10 dots, closely spaced, can be observed.

The situation will depend on how to want to treat the data. If you treat it as an image, putting white values (or whatever) between dots, you are already doing interpolation. In that case, a true 2D DFT would take a separable 1D DFT on both rows AND columns. So, you definitely need to venture into 2D, unless you will only  see horizontal frequencies; look at the last line of your image: there are 8 or 9 periods only, you are missing out the small slope.

I   cannot propose data interpolation without knowledge of the underlying process, and the $z$ value. My take is that you ought to try with true 2D non-uniform DFT first. Interpolation causes smearing.

In recent question (what is $1/f$ property of the visual world?), Frequency Synthesis of Landscapes was evoked, and a parametric spectrum modeling (eg decay in $1/f^p$) can provide you with parameters to compare different surfaces. For some spectral shapes, this can be naturally done in a Fourier domain. There should be related techniques under the field of "Surface roughness analysis".

My basic thinking for a procedure is the following: your different scans of surfaces may have different $(x,y)$ locations, and sizes. So it is not easy to compare their properties directly in the space domain, unless you regrid, interpolate, etc. which can cause harm. Non regular Fourier transforms can regularize the outcomes: from any two scattered grids, you can cast Fourier coefficients to the same 2D regular frequency grid, and infer global parameters from them: peaks, moments, parametric fits, etc.

So you can check NUFFT (NFFT, USFFT) Software and the related keywords:

  algorithms have been developed to overcome this limitation - generally
  referred to as non-uniform FFTs (NUFFT), non-equispaced FFTs (NFFT) or
  unequally-spaced FFTs (USFFT)

with a related paper: Accelerating the Nonuniform Fast Fourier Transform, by Leslie Greengard and June-Yub Lee. Apparently there are Matlab versions too:
NUFFT, NFFT, USFFT.

If this is not enough, you can try (lifting) wavelet transforms on 2D irregular grids, to catch the proper scale and local details. You could mesh it with MeshLab. And then use lifting wavelets for a multiresolution vision.


Spreading sequence and equalization based on a paper I am having difficulty in understanding how the Authors in this paper, An EM based method for semi blind identification of linear systems driven by Chaotic signals

have used the expressions derived from applying an EM_EKF estimation method when the driving process is a non-linear input. They have considered the noise model as AWGN. They then apply their method for in equalization. They apply the estimators derived in Eq(12)--Eq(15) to estimate the MA process with Rayleigh fading coefficients

The channel is a linear filter with taps whose magnitudes are modeled as Rayleigh random variables. Assuming, the input $(\mathbf{z}_n)$ take in 8 symbols from QAM constellation instead of two symbols $\pm1$ as in BPSK. The model is

$$ y_n = \mathbf{H} (\mathbf{z}_n) + w_n \tag{1} $$ expanding
$$y_n = h_1 z_{n}+ h_2 z_{n-1} + h_3 z_{n-2} + w_n \tag{2}$$

with channel coefficients as $\mathbf{h} = [h_1,h_2,h_3]$.

$$\mathbf{z}_n = [z_n,...,z_{n-p+1}] = [z_{n},z_{n-1},z_{n-2}] \tag{3}$$

$\mathbf{y}= [y_1,y_2,y_3, ...]$ is a one dimensional array (scalar valued) noisy output of n = 1,....,N data points.

The number of multipaths = 3. The Authors say that they have used a spreading sequnce in equalization. If so, would (1) still be the model and the expression for the estimators derived in the paper would not change? It is not clear form the paper how they have applied their method in equalization.

I am not aware what would be the functional form of the model in equalization application when a spreading sequence is used. Can somebody please explain how the method is applied in equalization using spreading sequence, the model (FIR) equation and how to apply the estimators using the spreading sequence. (An explanation with implementation would also be useful, if possible).

Thank you.

This is really a bit of a guess, since I haven't read the paper, but:

A spreading sequence can be designed/used to whiten a signal – that is the PSD of $y_n$ is the sum of a white process $w_n$ and a colored process $\sum\limits_{i=1}^3 h_i z_i$.

Based on observation of the receive signal you can estimate the PSD (e.g. by $\left|\mathrm{FFT}\right|^2$) and subtract your noise power density estimate (that you could get either from observing more bandwidth than your transmission, or by observing your band while you're not transmitting) for $w$; that way, you'd get the "spectral shape" of $h$, and thus, at least a start for equalization; for simple signals (like BPSK) equalizing the amplitude alone might suffice.

Now, for QAM etc., this might not really be the way to go.

However, we know the following about the autocorrelation of both your $w$ and your $z$:

Their autocorrelation is zero but for a single dirac impulse with the power of the signal at $\tau=0$.
Note that autocorrelation is a property of a stochastic process that you, again, cannot observe directly, but estimate from a sample autocorrelation. That's why when you look at graphs that correlate a signal with itself, you don't see that perfect behaviour but non-zero: your observation is of finite length.

Now, with that knowledge, a discrete-time non-multipath reception of a signal should have exactly the same property: with increased duration of observation, the autocorrelation of the receive vector should converge to a dirac impulse.

With multiple paths, the "shortest" path's receive signal should, however, also correlate with the the other paths – and thus, you should get "peaks" for every $\tau$ that is a path length difference.

That's pretty nice!

The "boring" way of finding those peaks is actually computing the sample autocorrelation – but since correlation has time complexity $\mathcal O(N^2)$, and you might even want to interpolate your $N$ samples first to get higher resolution, you'd do it the fast correlation with the FFT. Notice that you can then comfortably increase the resolution by zero-padding.

The algorithm for that would be
$$\begin{align}
 \mathbb R^K  \ni &\quad \mathbf{\tau}_\text{peaks} \\
 &=\text{largest}_K \big(R_{yy}(\tau)\big)\\
 &=\text{largest}_K \big(\mathrm{IFFT}_N\big\{\mathrm{FFT}_N\left\{y\right\}(f)\cdot\mathrm{FFT}_N^*\left\{y\right\}(f)\big\}(\tau)\big)\\
&=\text{largest}_K \big(\mathrm{IFFT}_{\text{hermitian},N}\left\{\left|\mathrm{FFT}_N\left\{y\right\}(f)\right|^2\right\}(\tau)\big)\\
\end{align}$$

The "$\text{hermitian}$" is but a hint – the FFT for real input data is faster to calculate in most cases than the general complex input case.

Notice how $\text{largest}_K \big(\mathrm{IFFT}\big)$ is a frequency estimation problem, where you're looking for $N$ frequencies – and looking for known-$N$ dominant frequencies especially well-solved by parametric estimators, where you can get very high resolution – without needing to pad the FFT in the first place. I'd like to point you to ESPRIT:
$$\begin{align}
 \mathbb R^N  \ni &\quad \mathbf{\tau}_\text{peaks} \\
&= \text{ESPRIT}_{N,K}\left(\left|\mathrm{FFT}_N\left\{y\right\}(f)\right|^2\right)
\end{align}$$


Real and Imaginary data from Network Analyzer I’m having a hard time grasping what exactly ARE the real and imaginary numbers in my s-parameter data. Most of the time I just set the NA up in logmag/phase to find the loss though, say, a filter or cable. My understanding is you can pretty much derive anything from the real/imag, especially s11. But are these numbers kinda useless (for lack of better term) until converted into other units?

The S paramater data are complex numbers that have magnitude and phase, or stated another way, real and imaginary components. The way we typically see the data presented on a network analyzer is with a Magnitude vs Frequency plot, and a Phase vs Frequency plot. If we capture the data we may see it listed as Real and Imaginary. The two are related as follows:

Given: Magnitude (K) and Phase ($\phi$)

$Sxx(\omega) = K(\omega)e^{j\phi(\omega)}$

(showing magnitude K vs frequency together with phase $\omega$ vs frequency since K at angle $\phi$ is $Ke^{j\phi}$ in compact exponential form.)

Has real and imaginary components that are also functions of frequency as follows:

$K(\omega)e^{j\phi(\omega)} = I(\omega) + jQ(\omega)$

Where I is often used for the real component (In- phase), and Q is often used for the imaginary component (Quadrature phase).

The relationship between the two is given by Euler's identity:

$K(\omega)e^{j\phi(\omega)} = K(\omega)cos(\theta(\omega)) + jK(\omega)sin(\theta(\omega))$

The real and imaginary numbers define the magnitude and phase versus frequency for the relevant S parameters. For example S11 defines the reflected signal with values between 0 and 1 (for passive devices) and phases anywhere from 0 to 360 degrees. An open circuit for example will reflect all incident power back toward the source in phase (so would be 1 (0 dB) at angle 0 degrees), while a short will reflect all the incident power back toward the source out of phase, 180 degrees (so would also be 1 but at angle 180 degrees).  S21 is the transmission from port 1 to port 2 of a 2 port network, and a simple cable is a good demonstration of phase: a lossless cable is just a delay, and a delay has a linear phase with frequency so would be magnitude 1 (0 dB) with a phase that increases linearly negative as frequency increases.

Thus to demonstrate the significance of the phase, you could predict the magnitude response of a comb filter that is created with two power dividers and a long length of cable between one leg and short on the other (split, go through the cables and then combine for a resulting two post network), by measuring S21 phase through each of the cables.

At certain frequencies the phase between the cables will be 180 degrees and thus S21 magnitude through the entire network will be nulled when the two signals combine, while other frequencies they will be in phase and the signal will be combined, thus resulting in a comb pattern in magnitude such as that shown in the normalized frequency plot below.

The phase of the reflected signal (S11) and transmitted signal (S21) is used often in many aspects of RF design including impedance matching and filtering design (since like in digital filter design filters can be constructed through the summation of weighted delays, and hence weighted phase shifts). Having a deep understanding of this is fundamental to RF design.


Issue with demonstrating amplitude modulation of sinusoid using MatLab I am trying to demonstrate Conventional and DSB-SC modulation of a signal $ \sin(2\pi1000t) $ using MatLab. My carrier signal is $ \sin(2\pi f_ct)$ where $f_c = 100\text{ MHz}$. Now according to theory, after modulation, the modulated signal should have a spectrum same as that of message signal shifted to the frequency of carrier on both sides. But I am getting the plot as that of the carrier signal in the frequency domain.
This is the code that I am using for DSB-SC modulation
% CODE FOR PLOTTING MESSAGE SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs = 20*100*10^6/3;
Fm = 1000; %Hz
t = 0:1/Fs:0.005-1/Fs;

Am = 1;

m = Am*sin(2*pi*Fm*t);

temp = size(m);
N = temp(2);

f = -Fs/2:Fs/(N-1):Fs/2;

M = fftshift(fft(m));
figure(1)
plot(f,abs(M))
grid('on')
title('Message signal')
xlabel('Frequency in Hertz')
ylabel('|M(F)|')

% CODE FOR PLOTTING CARRIER SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fc=100*10^6;%100 MHz

Fs = 20*Fc/3;
Ac = 1;%amplitude
c = Ac*sin(2*pi*Fc*t);

temp = size(c);
N = temp(2);

f = -Fs/2:Fs/(N-1):Fs/2;

C = fftshift(fft(c));
figure(2)
plot(f,abs(C))
grid('on')
title('Carrier signal')
xlabel('Frequency in Hertz')
ylabel('|C(F)|')

% CODE FOR PLOTTING DSB-SC MODUALTED SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s = Ac*m.*sin(2*pi*Fc*t);

temp = size(s);
N = temp(2);

f = -Fs/2:Fs/(N-1):Fs/2;

S = fftshift(fft(s));
figure(3)
plot(f,abs(S))
grid('on')
title('DSB-SB Modulated signal')
xlabel('Frequency in Hertz')
ylabel('|S(F)|')

And this is the spectrum of $S(f)$ in code. As you can see. there are only two peaks, one at $+f_c$ and one at $-f_c$(not shown in figure). But in reality, it should have been 2 peaks at $f_c+f_m , f_c-f_m$ and two peaks at $-f_c-f_m,-f_c+f_m$. So what is going wrong here? How can I fix this?

The problem is lack of frequency resolution in your plot. What you are seeing is in fact the two sidebands. Just zoom in to see them:

My advice is to use a much lower carrier frequency, such as 10 kHz. Doing so will make it much easier to study this modulation in Matlab, without sacrificing anything essential.


Converting speech audio to telephone audio I am trying to train a machine learning algorithm on telephony speech audio. However, there isn't really enough data for this anywhere that I can find. My solution is to just use speech audio from Youtube videos. How can I convert audio from a Youtube video to audio that simulates telephony audio?

Band-pass filtering with cut-off frequencies of 300 Hz and 3400 Hz should result in a good approximation. Try with a Chebychev filter or order not more than 6.

Then you may need to downsample your audio to 8000 samples per second, which is the standard for telephony.

P.S. The actual cut-off frequencies (especially the 3400 Hz) may be different according to different sources. I'll edit this answer when I find the official ITU recommendation.