texts: Faster way of getting 2D frequency amplitudes than DFT?

texts:
I'm making a small program which gets the DFT of an image to get a general idea of the image's overall orientation.  It does this by rotating a line in a radar sweep type pattern, keeping track of which angle has the highest sum of frequency amplitudes. (It actually only does $0-180^\circ$ though because the testing line is symmetrical across the origin)

Anyways, I was wondering, is there a less computationally expensive way to get frequency amplitudes of a 2D image than doing a 2D DFT?

I know that it is separable so can do better than the naive $\mathcal O(N^4)$ operation - I think it becomes $\mathcal O(N^3)$.  I also know that the FFT can do a 1D DFT in $\mathcal O(N \log N)$ instead of $\mathcal O(N^2)$.  Lastly, I know that it's massively parallelizable so could be (massively) multithreaded on either CPU or GPU.

Those things aside, are there any less computationally expensive ways to get frequency amplitudes of a 2D image than to do a full 2D DFT?


texts: Classification test by embedding true positives in high frequency recordings

texts:
I'm evaluating a cyclostationary approach to signal classification. My first step in doing so was generating signals and embedding them into a few specified parameter sets of the Watterson channel model. This seemed to be a good fit at the time. A colleague of mine proposed a new idea to me:

Description:

Embed/Add true positive digital signals generated by myself into broadband files recorded from an antenna?! Since I saw reason in his notion, I did as told. Using standard FFT-based signal processing I embedded digitally modulated signals into the broadband data files.

I subsequently found my true positive rate to be significantly less favorable compared to the Watterson channels tests. So his proposal, to me seemed to gap the difference from simulation to experimental study on a real HF channel. Since I now want to present my results to a scientific audience, I started musing over the influence of using real antenna data, with artificially embedded signals, as a test set:

The true positive signals on the one hand, acquired additive noise by embedding them into the real scenarios. But on the other hand, by not transmitting the signals, there is no significant influence resulting from  fading, multipath or spread distortion to my scenarios.

The cyclostationary approach predominantely requires a conservation  of the structure of the time signal. Therefore I assume, the relative strength of the detection didn't change through the additive noise?
The much more lively real data HF channel just introduces a lot of distorted similar signals which might be mistaken to be true positive signals by my algorithm?
Question:

How would you think does my choice of using real data with added true positives compares to real scenarios with truely transmitted true positives?


texts: how to calculate the maximum obtainable directivity for a microphone array?

texts:
Introduction

The relationship between the maximum obtainable directivity and aperture of an antenna is described by:

$$ D = \frac{4\pi A_e}{\lambda^2} $$

Where $D$ is the directivity, $A_e$ the effective aperture and $\lambda$ the used wavelength. The derivation can be found in the following document: pdf.

However, if the antenna does not scan as a sphere but in $x, y$ but not $z$. Following the same path as the above-mentioned document the equation becomes:

$$ D = \frac{2\pi A_e}{\lambda^2}  $$

The aperture of a lossless isotropic antenna is then also $A_e =  \frac{\lambda²}{2 \pi}$ instead of $\frac{\lambda²}{4 \pi}$

Problem

I cannot manage to use that equation for a microphone array using delay-and-sum beamforming.

For example, if I have an array of 25 microphones spaced as a grid and that my setup has the following properties:

Smallest distance $d$ between them is $0.01$ meters.
Amount of orientations $o$ that it looks at is $64$.
Frequency $f$ to be detected is 8 kHz.

The wavelength $\lambda = \frac{v}{f} = \frac{343}{8000} = 0.0429 $

The aperture $A$ is $A = \frac{\lambda^2}{precision} = \frac{0.0429^2}{\frac{2\pi}{64}} = 0.0187 $

Thus:

$$ D = \frac{2\pi A_e}{\lambda^2} =  \frac{2\pi A}{\lambda^2} = \frac{2 \pi \cdot 0.0184}{0.0429^2} = 62 $$

I rounded the numbers for the simplicity of formulating my problem, if you do not round then you get $64$ instead of $62$ which is the number of orientations that it looks at furthermore the results cannot be true because:

It does not involve the number of microphones used.
It does not involve the smallest spacing between microphones.
I probably made a mistake somewhere.

How can I calculate the maximum obtainable directivity for an array of sensors?


texts: Unexpected peak after up & down sampling

texts:
I am working on signals to check how approximation can change the results. I am testing on example that can change frequency from 22kHz to 20 KHz. I prepare this script in MatLab

clc;
close all;
clear all;
Fs1 = 22000;
Fs2 = 20000;
TimeStep = 1/Fs1;
Length = 100000;
time = (0:Length-1)*TimeStep;
InputSignal = sin(2*pi*8000*time);
nGcd = gcd(Fs1, Fs2);
wstop = Fs2/(Fs2/nGcd*Fs1);
wpass = 0.9 * wstop;
A = [1 1 0 0];
F = [0 wpass wstop 1];
FilterOrder = 255;
firCoeff = fir2(FilterOrder,F,A);
FilterOrder = length(firCoeff);
BigStep1 =  1000000000/Fs1; %this is 45455,55
BigStep2 = 45455;
Buffer24 = Timestamps(Length, BigStep1); %buffers of timestamps
Buffer24wrong = Timestamps(Length, BigStep2);
SmallerStep = 50000;
Length2 = Length*BigStep2/SmallerStep;
Buffer20 = Timestamps(round(Length2), SmallerStep); %output timestamps
upY = upsample(InputSignal,Fs2/nGcd);
ApplySignal = applyFilter2(FilterOrder,firCoeff,upY,Buffer20,Buffer24);
FilteredSignal= downsample(ApplySignal,Fs1/nGcd);
Length = length(FilteredSignal);
NFFT = 2^nextpow2(Length);
Y = fft(FilteredSignal,NFFT)/Length;
f = Fs2/2*linspace(0,1,NFFT/2+1);
FilteredSignal2 = applyFilter2(FilterOrder,firCoeff,InputSignal,Buffer20,Buffer24wrong);
LengthW = length(FilteredSignal2);
NFFT2 = 2^nextpow2(LengthW);
Y2 = fft(FilteredSignal2,NFFT2)/Length;
f2 = Fs2/2*linspace(0,1,NFFT2/2+1);
figure
subplot(2,2,1)
plot(f,2*abs(Y(1:NFFT/2+1)))
hold on
title('Wrong Input with up&dw sampling')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
subplot(2,2,2)
plot(f2,2*abs(Y2(1:NFFT2/2+1)))
hold on
title('Good inpuit olny without anysampling')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
subplot(2,2,[3,4]);
Index1 = GetIndex(Buffer24,Buffer20(5530),1);
Index2 = GetIndex(Buffer24,Buffer20(5540),1);
plot(Buffer24wrong(Index1:Index2),InputSignal(Index1:Index2),'-o',Buffer24(Index1:Index2),InputSignal(Index1:Index2),'-o', Buffer20(5530:5540),FilteredSignal(5530:5540),'d');
hold on
set(gca,'XTick',Buffer20(5530):SmallerStep:Buffer20(5540))
grid on
legend('Wrong signal','Up/Dwsampled Signal', 'Output Timestamps')
title('Input Signals')
xlabel('Time')
ylabel('|Y(f)|')

An also function GetIndex:

function Index = GetIndex(EquidistancyBuffer, ReqTime, StartIndex)
LastIndex = length(EquidistancyBuffer);

while((LastIndex-StartIndex)>1)
    Mid = round((StartIndex+LastIndex)/2);
    ActualTime = EquidistancyBuffer(Mid);
    if (ReqTime==ActualTime)
       LastIndex = StartIndex;
    elseif (ReqTime<ActualTime)
        LastIndex = Mid;
    else
        StartIndex = Mid;
    end

end
Index = Mid;

end

Function applyFilter2:

function yOut = applyFilter2(FilterOrder, Coeff, Signal, ReqTimestmps,EquidistancyBuffer)

StopIndex=length(ReqTimestmps);
for i = FilterOrder+1: StopIndex
    if(i == FilterOrder+1)
        startIndex = 1;
    else
        startIndex = i-1;
    end
    Index = GetIndex(EquidistancyBuffer,ReqTimestmps(i),startIndex);
    yOut(i-FilterOrder) = OutputCalculation(Index-255,FilterOrder,Coeff,Signal);

end

and Timestamps:

function Timestamps_Buffer = Timestamps(Size, Step)

for i = 1: Size
    Timestamps_Buffer(i)=(i-1)*Step;
end
end

I expected on the left subplot only peak on almost 8kHz but I have one closer to 9KHz. When I change Bigstep2 to equal to BigStep1 the results are almost the same. Please explain me what I did wrong. I attached the result of calcultaion also.


texts:
I have a rtlsdr dongle which apparently has a max sample rate of about 2.4 MHz.  I'm wondering how it's possible to capture signals at say 900 MHz with this device.  Doesn't the sample rate need to be 2x900 or 1.8 GHz to capture without aliasing?

Sorry, I'm so new to this stuff!

Thanks!


texts: Instantaneous frequency vs fourier frequency

texts: Lets consider a pure sine signal at $\nu$ that is chopped using square pulses (like a burst mode on signal generators). My understanding is that instantaneous frequency is $\nu$ when oscillations are ON and 0 when they are OFF. On the other hand fourier spectrum is constant over time and contains also other frequencies, since it is not pure sine anymore. Is this correct? which one is used when calculating some frequency dependent physical quantity?


texts: How to find frequency response, stability, and causality of a linear system?

texts:
I have the following transfer function:

$$H(s)=\frac{s}{(s+1)(s+2)}$$

How can I find the gain and phase response of the above system? I know the first step has something to do with substituting $s = j\omega$ into $H(s)$. How can I find whether the system is stable and whether it is causal?


texts:
I am trying use the FFT in a different way then most people ask about. I want to be able to take a picture of a graph with regular repeating vertical lines, and to process the image to determine how far in pixels the lines are apart on average. I have tried canny edge detection and hough line detection and I don't think I can optimize the images enough to accurately detect only the lines I am interested in.

So, my attempt is to scan 10 lines of the picture and to accumulate the pixel values into bins corresponding to the pixel column. What results when you graph it is a very nicely appearing waveform. When I perform a DFT or FFT on this, I can find a peak that I believe should be the frequency of the line repetition. (This may be a faulty assumption)

My question is, what does this number correspond to? i.e. I think I am confused with what my sampling rate would be because it is in pixels. I do think that this is a valid use of the FFT, but am falling right here at the point when i think I should be successful.

As an example. I created a picture that is 300 pixels in width. There are 1 pixel width lines drawn at exactly 30 pixel intervals. I found 2 peaks, one at 75 and one at 225 (which seem symmetrical) for the real component. (I do not believe that the imaginary component should play in??) I know the lines are 30 pixels apart. How does the 75 and or 225 relate?

I am trying really hard to get this, and I am grateful for any help you could recommend. At this point, I am giving up on edge detection, and want to try this approach.

Thank you in advance.


texts:
I am attempting to implement a basic bandpass filter but the center frequency of my filter seems to be irrelevant as I can change it all I want but it has no effect on the filter.
where am I going wrong?

I have given my calculations here in the hope that someone will be able to point out where I'm making a mistake.

Filter specs:
supression of DC and Fs/2 (zeroes at +1 and -1)
center frequency at pi/4
bandwith of pi/16

2(1-R)=pi/2
1-R=pi/32
R=(pi/32)-1
R=-0.9018252296

 K*(z-1)*(z+1)
---------------
z-R*cos(pi/4)+R

     K*(z^2 -1)
-------------------
z^2-R^2*cos(pi/4)+R

all leading to a difference equation of

Y[n]=R^2*cos(pi/4)*Y[n-1]-R*Y[n-2]+K*(X[n]-X[n-2])


texts:
I have two questions regarding systems with no poles:

Why does a system with no poles have a finite support?
Why if the system has a finite support it means that it is BIBO
stable?


texts: The magnitude spectrum of a sharpening filter

texts:
I'm trying to derive an expression for the magnitude spectrum of the following sharpening filter.
$$
g(m,n) = \delta(m,n)+\lambda (\delta(m,n) - h(m,n))
$$
where $\lambda$ is some positive constant and $h(m,n)$ is a low pass filter.
$$h(m,n) = \begin{cases}
\frac{1}{25}, & \text{for $|m| \leq 2, |n| \leq 2$} \\
 0, & \text{otherwise}
\end{cases}$$
Using the linearity of Fourier transform, I get
$$ G(u,v) = 1+\lambda(1-H(u,v))$$
where
$$
H(u,v) = \frac{1}{25} \cdot \frac{\sin\frac{5\mu}{2}\sin\frac{5\nu}{2}}{\sin\frac{\mu}{2}\sin\frac{\nu}{2}} \cdot e^{j\cdot 6(\mu+\nu)}
$$
So, to get the magnitude spectrum, I make the following steps
\begin{align*}
|G(u,v)|^2 &= \big[(1+\lambda)-\lambda H(u,v)\big] \big[(1+\lambda)-\lambda H^*(u,v)\big] \\
& = (1+\lambda)^2-(1+\lambda)\lambda\cdot 2\Re \big\{H(u,v)\big\}+\lambda^2|H(u,v)|^2
\end{align*}
where
$$
\Re\big\{H(u,v)\big\} = \frac{1}{25} \cdot \frac{\sin\frac{5\mu}{2}\sin\frac{5\nu}{2}}{\sin\frac{\mu}{2}\sin\frac{\nu}{2}} \cdot \cos(6(\mu+\nu))
$$
Therefore, to obtain the magnitude spectrum I should take the square root of the expression above.
However, the picture I get when I'm trying to plot it does not coincide with something that I'm expecting to see. I mean it's a high pass filter, hence there should be something like that.

But I get this.

Maybe there is a way to simplify this expression or maybe there is something that I'm doing wrong.
Thank your for your help!


texts: What is hidden in DFT below frequency resolution?

texts:
Setting: signal sampling rate 1000Hz, DFT 1s 1000 samples. This yields frequency resolution 1Hz.
I know, this is the limit of meaningful frequency and when I want to read 0.5Hz then either sampling rate should be doubled or samples acquisition doubled.
What bothers me is that when I compute DFT of 0.5Hz in given setting on artificial signal (sparse signal, some k Hz cosines and some 0.5 Hz cosines, e.g. 2Hz, 5.5Hz, 12Hz, 14.5Hz etc.) without any spectral window, the results seems valid.
I know this shouldn't happen, with dense signal, non-perfect frequencies (multiples of frequency resolution) and noise situation gets complicated.
What exactly I read below frequency resolution and what are conditions it may work?


texts: Radon Transform reducing the number of parallel beams (MATLAB)

texts:
I am trying to understand different filters' and other things' effects on Radon Transform by using MATLAB. First I upload an image, then take its Radon transform with radon function and then reconstruct it with iradon function.

Now I wanted to see the difference between undersampling the projection data by reducing the number of parallel beams by a factor of 4 and reducing the number of rotation angles by a factor of 4. However according to what I understood from radon and iradon functions' descriptions, I can only change the number of rotation angles. Is not it true? If not, what is the way to do it?

Here is what I do to decrease the number of rotation angles by a factor of 4:

I = imread ('asdasd.tif');
Rad = radon (img,0:4:179);
I_ram = iradon(Rad,0:4:179,'linear','Ram-Lak');
I_ham = iradon(Rad,0:4:179,'linear','Hamming');


texts: Deconvolve a symmetric sequence into minimum-phase and maximum-phase sequences

texts:
I learned in school that any real-valued, even-symmetric sequence can be expressed as the convolution of a minimum-phase sequence and its time-reverse, which is maximum-phase.  In my case, the real-valued, even-symmetric sequence is $sinc(x)$.  What causal, real-valued, minimum-phase sequence convolved with its time-reverse produces $sinc(x)$?

I am familiar with the traditional approach of computing the frequency-domain zeros and constructing the minimum-phase sequence from those inside the unit circle and half of the pairs on the unit circle.  However, in this case, the sequence has infinite length.  I would prefer a closed-form expression for the minimum-phase sequence function.

In practice, I will be using a windowed version of the sequence -- an FIR filter.  I would also like to learn how to deconvolve any real-valued, even-symmetric, finite sequence (FIR filter) into its minimum-phase component, whose frequency response magnitude is the square-root of the symmetric FIR.  I would like a time-domain method for doing this because there are precision problems with computing frequency domain zeros in very long sequences, such as a million taps.


texts: Creating 2D array to play audio in left and right channel in matlab

texts:
I'm new to matlab sound processing and I am trying to create a 2D array for the sound code below. I have to make the sound have 400 samples(i.e 10ms) at 20hz, but having trouble getting the code to work. Any help would be greatly appreciated.
The steps are:
(1)    Construct a two-dimensional array (n x 2  array, i.e., 2 columns, n rows, n is the number of samples per sound channel).  The array is the sound signal, the first column is the sound for the left channel (left speaker), the second column is the sound for the right channel (right speaker)
(2)    Then play the sound.

Fs = 44100;
t = [0:4/Fs:4-4/Fs];
freq = 20;
f1 = sin(2*pi*freq*t);
f2 = sin(2*pi*2*freq*t);
%Sound 1
sound(f1,Fs);
%Sound 2
sound(f2,Fs);
%Play Consecutive
f12 = [f1 f2];
sound(f12,Fs);
%Play together
f_12 = [f1+f2];
sound(f_12,Fs);


texts: Is sound the best way to transmit a signal under water?

texts:
Could sound be used in transmitting sound waves under water digitally? Could the sound then be converted by a repeater to transit by radio waves? What other ways can a signal be transmitted other than sound over long distances under water?

https://space.stackexchange.com/questions/32421/how-would-a-transmission-be-relayed-through-bottom-of-europa-ocean-to-earth


texts: Generating LDPC codes of arbitrary sizes in MATLAB

texts:
I am having problems generating LDPC codes of arbitrary sizes in MATLAB primarily because MATLAB insists on the parity check matrix that is input into its ldpcenc function has the sub-matrix that encompasses the last $n-k$ columns being invertible. The error that is thrown by the ldpcenc function when an unsuitable matrix is input is:

The last (N-K) columns of the parity-check matrix must be invertible in GF(2)

I am generating parity matrices with the online tool found HERE which accepts the maximum node degrees and the code rate as input and outputs optimised parity check matrices for BEC and AWGN channels. An obvious solution to this problem would seem to just convert my parity check matrices to standard systematic form - I tried doing this by using THIS function from Matworks file exchange but the resulting parity check matrix (even though it gets accepted just fine by the ldpcenc function) doesn't seem to 'work' (i.e a code with rate = 1/2 and block size 20 can't correct the received code even if there is just 1 bit in error).

How can I generate LDPC codes of arbitrary sizes in MATLAB? Any help appreciated.


texts:
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.


texts:
I'm trying to reconstruct an audio signal, and below is what i got after rescaling.

How can i obtain centered audio signal such as below?

SOLVED:

ruoho ruotsi's answer addresses my problem :)

i found someone asked the similar question before

How to remove or filter the drift problem in measured Strain signal?


texts: Choosing clock frequency and number of bits of an OFDM transceiver

texts: I want to implement an OFDM transceiver chip and my question is how to choose the clock frequency of the IFFT block? The number of output bits of the IFFT block? Also, how has the clock frequency and the number of output bits of the digital-to-analog converter (DAC) to be chosen?


texts:
I would like the chirp to end in phase zero.
Chirp time or end frequency may vary slightly.
Now I'm checking the maplot output.
I observe the same thing at the output of the sound card.
In this solution, the chirp does not always end in phase zero.
help find a solution.
Thanks for the help.
Andrew
from pylab import *
import numpy as np
from scipy.signal import chirp

f0 = 7000
f1 = 17000
samplerate = 192000
T = 0.0013
T = np.ceil(T*f1)/f1 # new T
t = np.arange(0, int(T*samplerate)) / samplerate
w = chirp(t, f0=f0, f1=f1, t1=T,phi=270, method='linear')

fig, ax = subplots(figsize=(6,1))
ax.set_title("Chirp ")
ax.plot(w)
show()


texts:
Q1. Consider a system implementing a rational sampling rate change by 5/7: for this, we cascade upsampler by 5, a lowpass filter with cutoff frequency π/7 and a downsampler by 7. The lowpass filter is a 99-tap FIR.
Assume that the input works at a rate of 1000 samples per second. What is the number of multiplications per second required by the system? Assume that multiplications by zero do not count and round the number of operations to the nearest integer.
Please tell the approach on how to start this question.
The generic order of sampling rate change is "upsample", "filter", "downsample". That means that we need to run the filter at the upsampled rate, i.e 5000 but given that it is 99tap filter so shall i directly multiply 99*5000=495000. is it okay wrt the asked question or am i missing something


texts:
Given something like 103 data points ($N=103$), a DFT will give back 103 frequency values. Then to do something like filtering the high frequencies involves setting the high frequency values from the DFT to zero, and doing the inverse DFT to get back 103 data points that represent the original signal without the high frequencies.

This makes sense to me given all 103 data points all at once. But what about streaming an audio WAV file that is rather large (say for example $N=10^5$). If one wanted to filter the high frequencies then the approach i just described on the entire data of $10^5$ points logically makes sense. But that is not reasonable when streaming a WAV file for playback. What is done to high frequency filter a streaming playback of an audio file?


texts: comparison between Fourier transform and autocorrelation function

texts: If the Fourier Transform of a signal is enough to determine its spectral density, why do we complicate our lives by defining the autocorrelation function?


texts: Can't understand multi intensity valued mask

texts:
I am trying to do nuclei detection from a histopathological image. The original image and the mask (ground truth image or segmented image) are as follows:

The tissue image is (128,128,3) [RGB] and the mask is (128,128). Now, if I visualize with grayscale colormap:

My ultimate aim is to get bounding boxes for each nuclei in the mask.
What I have understood till this point: Do segmentation on tissue image, obtain the mask, find the bounding boxes on mask. Then I questioned why to get bounding boxes from mask, why not the original image? Then after some searching it led to conclusion that masks are more refined images only emphasizing required portions and setting remaining to 0.
Now, from my understanding a mask would be a binary image with certain pixels as 255 and rest as 0, so that when it is bitwise AND-ed with the original image we get only the required portions (wherever pixel values were 255).
So, in the given mask, how come there are different intensity values for nuclei? Is it due to the fact that some nuclei are strongly represented and some are weakly represented?

If so, then if I threshold the mask to 255 for every pixel value greater than 0, it should ideally become the mask which I described in my understanding, right?
If not so, then what is the meaning of this multi-intensity valued mask? It can't be applied like a binary mask for AND operations (since it has values like 38,60,110 apart from 255).
Even though it's a one channel image, the default colormap sets nuclei to different colors in the mask. Does it indicate that the mask is a m-ary mask with same color nuclei belonging to a certain class/label and if AND-ed with original image, then would yield only those nuclei of that class?


texts:
Okay, so I am trying to address a biological problem here. I have asked a variant of this question here, but I was not able to get to solving it.

The question:

Four objects are vibrating constantly with some variable frequency and amplitude. I have the X-Y coordinates of the objects w.r.t time. Now, I want to find if there is any correlation between the objects movements.

Also, if there is any correlation between the movement of an object with itself w.r.t time.

I am not sure what to do, I have absolutely no idea about the mathematics involved, but can try to understand what needs to be done.

Thanks!
A


texts: How to combine a rotation matrix and a stretch matrix into a single matrix for easy Fourier Transform

texts:
For full disclosure, this is related to homework.  I have to find the Fourier Transform of a function that I've boiled down to the following.

I have a function $f(x,y)$ that I can think of as another function $g(x,y)$ plus a stretch and rotation.  (For example, think of turning a circle into an ellipse).  So, I have a stretch matrix of:

$$\pmatrix{ \frac{1}{A}& 0\\\ 0& \frac{1}{B}}$$

And I have the standard rotation matrix:

$$\pmatrix{ \cos \theta& -\sin \theta\\\ \sin \theta& \cos \theta}$$

I know that for computing the Fourier transform, $f(Ax)$ is the same as $\frac{1}{det A} F(A^{-T}u)$.  So, if I can combine the stretch and rotation into one matrix A, I'm home free.  How can I combine these two matrices?

Thanks!


texts: PSD - Correct Formula for removing Noise

texts:
Ok, so finally some things are making sense.

I'm trying to scale the results, which I have done by the following equation:

$$
PS_{rms} = \frac{2 \cdot \left |y _{m} \right |^{2}}{
   S_{1}^{2}
} m = 0.....N/2
$$

Where   $y_{m} = (re*re+im*im)$

This produces the following:

(Sorry about the axis - I have not had time to do these)

Prior to scaling the results, I got the following:

Essentially, so far, so good (I hope!)

Now I want to remove the noise around the signal, in theory, so that the main peaks are just showing and remove all the peaks that are not. So in theory, I would just be left between points (60 - 80).

I have been trying to read, and, understand the "Section 10: Averaging and overlapping" in this paper:
http://www.rssd.esa.int/SP/LISAPATHFINDER/docs/Data_Analysis/GH_FFT.pdf

It assumes that the usual remedy is to take the average of $M$ elements and reduce the standard deviation of the averaged result by a factor of $1 / \sqrt{M}$ however, what would be $M$ in this case? How could I therefore reduce the standard deviation, in order to remove noise from the result of scaling the PSD?

This is the only part to which I am now confused about, I essentially need to reduce the noise levels at the start of the result, and, at the end of the result thus giving me only the important peaks within the signal.


texts:
I have a .wav file of human speech, sampled at 48kHz, with useful signal up to 5kHz, and almost no noise. If I do Hilbert transform in MatLab, I get the Analytic signal:

$$ s_\mathrm{a}[n] = s_\mathrm{real}[n] + j \cdot s_\mathrm{hi}[n], $$

where $s_\mathrm{hi}[n]$ is a Hilbert transformation of $s_\mathrm{real}[n]$ from the .wav file. Now for each sample I can calculate the phase:

$$ \phi[n] = \operatorname{Arctan}\left( \frac{s_\mathrm{hi}[n]}{s_\mathrm{real}[n]} \right) $$

and take time derivative from it. That will be the speed of phase changing. I don't want to call it "instantaneous frequency", that's too controversial.

Is any smarter way of doing this in MatLab? Probably this:

Fs/(2*pi)*diff(unwrap(angle(hilbert(s_real))));

I did this in Adobe Audition plugin, in Delphi. Way too many spikes, every time when the envelope $\sqrt{s_\mathrm{real}^2[n] + s_\mathrm{hi}^2[n]} $ was around zero.


texts: How does the band pass wavelet act as a high pass filter in DWT

texts:
I am understanding that wavelets have a bandpass filter response due to their compact support, and when doing decomposition they are used as a high pass filter along with the scaling function as a low pass filter in DWT.
So my question is how does the wavelet perform as a high pass filter and does not attenuate near nyquist frequencies since it is bandpass in nature?


texts: Is it useful to apply edge detection filter as a preprocessing step before applying CNN?

texts:
Background: my future boss would like me to apply a number of low level filters to industrial images before applying machine learning methods on them.
From my knowledge this doesn't make any sense as CNN will learn to do exactly the same, making the process redundant, in the best case, and counterproductive in the average case, as it will alter the natural information, reducing the original entropy. Moreover this could introduce technical problems in the learning chain as the version of the learnt ANN will depend on the version of the preprocessor, etc.

Is there any reason to apply such kind of filters before CNN?


texts: High Pass Filter, effect of applying it on this image.

texts: I am having question about the HPF. In particular, it is having a high pass guassian filter which is applied to the original image below. High pass filter should only detect edge (lots of change in terms of intensity of the pixel). But how come the resulting image, still contains the two "white" things in the middle? Since they are not edge (no change in terms of intensity), should they be filtered out as well?


texts: I'm a Software Engineer Undergraduate. I'm trying to detect drum strokes(beats) in real-time in python. I used amplitude threshold-based onset detect (beat detect) algorithm but it takes other noises also onsets and it omits some actual beats with low amplitude. Can you have suggestions to optimize my code or another way to solve this?


texts: Reconstructing/interpolating small regions of a bandlimited signal by taking the fewest possible samples

texts:
I have a signal which is bandlimited and can be sampled at arbitrary continuous positions. The value at any position is given by an expensive computation. I need to do some further computation on various arbitrary small regions of this signal.

My plan is to take advantage of the known bandlimiting to sample the signal at or slightly above the Nyquist rate within these regions that I need to process, and then do my computation on a reconstructed approximation of the region. This reconstruction needs to be in a form that allows for evaluation at arbitrary time t within the reconstructed range. Because of the expense of computing each sample, I need to minimise the number of samples required to reconstruct each region. The regions are sized such that they are between 8 and 16 samples wide at the Nyquist rate.

The thing that makes this different from standard usages of sinc interpolation is the small number of samples. That's the focus of this question. It's typical for sinc based reconstructions to experience problems near the edges of the signals, but when a signal is only 8 samples wide, most of it is edge.

My question concerns how to optimally reconstruct these small regions with minimal samples. I particularly would like to know whether it's possible to achieve a reasonable reconstruction with all samples being within the region I'm reconstructing. I suspect that the answer is no, because there's insufficient information to reconstruct towards the edges of the region, but I would like to have this confirmed. My experiments with various windowed sinc methods support this suspicion.

Assuming that I can't rely on only samples within the region, is there anything better I can do than taking a few extra samples on either side of the region and then continuing to use basic sinc methods?

I don't have a firm definition of a good enough reconstruction, but for the purposes of the question, let's say that at all points along the reconstruction, the error should be no more than 10%.


texts:
I have a car with 8 lidars, each with a field of view of 60 degrees. My car looks like this

How can I merge all the lidar readings into 1 point cloud?


texts: How to create a swept / chirp signal between a given starting audio signal and a given ending audio signal

texts:
I can create a swept / chirp signal using matlab / octave with the code below (also see spectrum taken).  But how can I create a swept / chirp signal between a given starting audio signal and a given ending audio signal?  (One way to look at it is morphing one image into another image but using two audio signals)  or a sliding / chirp phase vocoder

clc,clear all
swept_startfreq = 500;
swept_endfreq = 225;
fs = 8000;
swept_dursec = 15;  % duration of signal in seconds
swept_t = (0:swept_dursec*fs)/fs;  % Time vector

if (swept_startfreq == swept_endfreq)
      swept_phase = swept_startfreq * swept_t;
else
    swept_alpha = abs(log(swept_endfreq/swept_startfreq))/swept_dursec;
    swept_b = (swept_startfreq - swept_endfreq)/(1 - exp(-swept_alpha*swept_dursec));
    swept_a = swept_startfreq - swept_b;

    swept_phase = swept_a*swept_t - (swept_b/swept_alpha)*exp(-swept_alpha*swept_t);
endif

swept_sig_out = .8* cos(2*pi*swept_phase); %equation used to create signal
wavwrite([swept_sig_out'] ,fs,32,strcat('/tmp/swept_sig_out.wav'));  % export file

The frequency spectrum of the original code / chirp or swept signal it creates

So I would have an initial audio signal (use as the initial audio signal for testing) and an ending audio signal (use as the ending audio signal for testing).  How can I do a "morphing" between these two signals?

Here's what an example of the final file would sound like example of final audio signal when processed (please not this isn't exact but it shows how the pitch gradually changes from the original signal to the end audio signal over time  If you would like to know how I created the final file I used the initial signal and Audacity with the effect - Sliding Time Scale/Pitch Shift to get as close as possible to what the ending signal sounds like.

Added from comments: Start with the frequency with the largest amplitude in the initial signal and "morph" it towards the frequency with the largest amplitude found in the ending signal

It would sound like "one one one one one" gradually decreasing or increase in pitch depending on the ending audio signals largest amplitude frequency found.   Please note the Initial and ending audio signal may not be the same type of signal, the initial signal could be "one one one one" and the ending signal could be "two two two two"

Ps:  I'm using Octave 4.0 which is like matlab


texts: How poles are related to frequency response

texts:
I have recently fallen into fallacy, considering pole s=1 as there is infinite response at frequency 1. Yet, response was only 1. Now, can you derive the frequency response, given the poles?

Secondly, the theory says that a system is stable when poles are in left s-plane and, thus, decay in time. But, wait. Does'n "pole" mean the infinite response -- the growth in time?

Finally, is it right question in DSP? IMO, D stands for digital whereas s-domain is analog. I do not find s-plane or Laplace transform tags to label my post.

update Thanks for the answers. It seems that I have got it except the one minor but fundamental thing -- the relationship of poles (and zeroes) with frequency. Basically, why are eigenvalues (or, how do you call the $s$ operator/variable) related with frequency? It should be somehow related with exponential growth and Laplace transform. I quite understand that poles happen to be eigenvalues (especially for discrete recurrences). But, how is this related with frequency?


texts: why human voice is band limited while it's time limited?

texts:
According to Wikipedia, time-limited signals can not be band-limited and vice versa. Also we know that human voice is band-limited and has frequencies from $50\textrm{ Hz}$ to $4\textrm{ kHz}$.

So how is the human voice band-limited while it's time-limited ? (consider a man speaking in a period of time to have a time-limited voice signal)


texts: How can I convert audio to a series of notes using pitch detection?

texts:
I have a good pitch detection system set up, and I would like to return a series of notes given an array of audio samples.

My current approach is as follows: I have a moving window across the audio signal, and I calculate the pitch of each window. Afterwards, I segment the audio into different notes by detecting the silent regions (i.e. where the pitch detector returns null). I then simply take the average of each note region.

Unfortunately, this has not been giving me such good results. The pitch detector does seem pretty accurate, but the issue is it doesn't seem to segment the notes very well. It only really works when I leave a long pause between each note when I record the audio. I would like some way for it to detect a change in notes without having to rely on a large silent region.

Any ideas would be greatly appreciated!


texts:
I have to identify the individual keys presses from a DTMF signal. There are three keys presses together as seen in the image. The signal has a sampling frequency of 8kHz. Each tone lasts between 0.1 and 0.2 of a second and there is a gap between tones of at least 0.05 of a second. However, there is random jitter in the timing so the actual timing will be unpredictable. The level of noise on the output is given by the signal to noise ratio of, SNR. I am unable to proceed with this.
Any help with this will be really helpful


texts: How do I show an impulse response leads to a zero-phase frequency response?

texts:
I'm trying to understand how to show that with real coefficients, the phase response of a filter is 0. Here is the impulse response

$h[n] = b_1d[n+1] + b_0d[n] + b_1d[n-1]$

How should I approach showing this?


texts: Filter to speed up convergence of geometric series

texts:
I want to speed up convergence of geometric series
$\ \  lim + c^n \ \ $ where $|c| < 1$ and $lim$ is not known.
A basic low-pass filter, here Butterworth, is good for $-1 < c < 0$ but not $0 < c < 1$:

butter  N 2  Wn 0.5:  [.293 .586 .293] / [1 0 .172]
|response|: [100 100  99  97  88  71  47  25  10   3   0]
(-.8) ^ n  [100 -80  64 -51  41 -33  26 -21  17 -13]  ->  [29 35 -4 -7  1  1  0  0  0  0]
.8 ^ n     [100  80  64  51  41  33  26  21  17  13]  ->  [29 82 90 62 45 38 31 25 20 16]

Is there a filter that does both ?

Why: many iterative processes that generate points $X_n$ in 1d, 3d, 100d
converge to an $X_{lim}$ slowly,
but roughly geometrically:
either one-sided like $.8^n$, or oscillating like $(-.8)^n$ .
It would be nice to have a simple, robust filter that speeds up both.

Shanks transformation
is a nonlinear transform that works really well on geometric series.
But on linear input 1 2 3 it'll divide by 0.
(Check for that, push on regardless ? practical, but unsatisfying.)
(I believe that
Aitken's delta-squared process
is the same as Shanks -- not sure.)


texts: How to calculate SNR variation with period?

texts:
I want to calculate the SNR variation of a signal with respect to the period. The current technique I am using is calculating SNR iteratively by applying narrow bandpass filters of different range, but the computational time for this is painstakingly high.
Is there a way to use fft or something like that to calculate the variation of SNR with period?
(I asked this question in stack overflow and I got a suggestion to ask here)


texts: Python equvialant of zero-phase IIR band-pass filter

texts:
A continuation on Matt L. answers, where

Zero phase filtering with IIR filters is achieved with forward-backward filtering, as implemented in Matlab's filtfilt function.

Main question: How can I apply zero-phase IIR band-pass filter on Python?
SciPy has scipy.signal.filtfilt function, but I am unsure on the configuration to perform a zero-phase IIR band-pass filter from 5 Hz to 50 Hz.


texts: turn circular convolution into linear convolution by zero padding: A special case

texts: We know that, multiplying a kernel and signal spectrum in Fourier domain will lead to a circular convolution and not a linear convolution, so in order to it become linear convolution we must zero pad your signal and kernel before taking the Fourier transform (up to M+N-1 where M is the signal's length and N is the kernel's length). However, I have already known the analytical solution to the Fourier transform of the kernel, while I don't know the analytical form of the kernel. I know that I can inverse discrete Fourier transform $\mathcal{F} \{kernel\}$, then sample and do the zero-padding, but I will lose the advantage of the "analytical form" in this case. Is there any good way for me to turn the circular convolution into linear convolution?


texts: What is the inverse Laplace transform of squared denominator term?

texts: Referring to the image below, what would the inverse Laplace transform be? I can't seem to find any tables that include this case.


texts: Heisenberg boxes image for Matching Pursuit

texts:
I'm using the MATLAB algotirhm Matching Pursuit to obtain the Gabor atoms forming a signal.
But now I want to graph them in a "Heisenberg boxes" graph.
Does someone how if is there a function to graph them like this?:

Or maybe I should write a function for graph. In that case, I know the time and frequency which are the coordinates of the atom, also I have the scale which give me the width in time. But how can I know how much wide is in frequency? I suppose that I must have a certain area and by knowning the width in time I can calculate the width in frequency but I'm not sure.

-------------------------UPDATED----------------------------

I made my own function:

function plotGaborAtomos(N,params,fs)
%N is the length of the atoms
%params gives the parameters from Gabor atom
%params is size 3xnCols where first row has the scale, second row has timeShifts and third row the freqs
%params = scales timeShifts freqs
nCols = size(params,2);
a = [ 0 N/fs 0 fs/2];
plot(a(1)-.5, a(3)-.5);
axis(a);
hold on;
titlestr = 'Sonido prueba';
titlestr = ['Heisenberg boxes: ' titlestr ];
title(titlestr);
xlabel('Time');
ylabel('Frequency');

area = N;
for i=1:nCols
    s = params(1,i);
    tCentral = params(2,i)/fs;
    fCentral = params(3,i)*fs;
    t0 = tCentral - (s/fs)/2;
    t1 = tCentral + (s/fs)/2;
    wideFreq = area/(s) * fs/2/N;
    f0 = fCentral-wideFreq;
    f1 = fCentral+wideFreq;
    plot([t0 t1 t1 t0 t0], [f0 f0 f1 f1 f0])
    plot( tCentral, fCentral, 'o')
end

Using it in a signal of length 256 gives me this:

As can be seen, the sampling frequency is 44100Hz


texts: gaussian window for designing filters in DSP

texts:
I understand that we do not use rectangular window as a window in filter design because of the ringing/gibbs artifacts due to the sidelobes of the resulting sinc function in the time domain.
Consequently, most DSP books mention hanning, hamming windows because they reduce the ringing artifacts.
But how about using a simple gaussian window? It is symmetric, invest DFT will also be gaussian which will have much lesser sidelobes than sinc function, etc.
Is the ringing artifact less for hanning window than the Gaussian window? For me a narrow transition band is not the topmost priority.


texts: IDTFT of convolution in the frequency domain

texts:
I have tried everything. If you actually know how to solve this could you provide a hint?
$$ e^{-2j\Omega}\frac{ \sin\left( \frac{7\Omega}{2}\right)}{ \sin\left( \frac{\Omega}{2} \right)}\star \frac{\sin\left( \frac{10 \Omega}{2} \right)}{\sin\left( \frac{\Omega}{2} \right) }$$
Ideally I would like to find the Fourier of every "fraction" in separate and then use properties:
$x(n - n_o) \rightarrow e^{-jn_0\Omega}X(\omega)$ so I don't mind for $$ e^{-2j\Omega}$$
but I have 2 problems:

I cannot use $\displaystyle \frac{\sin\left(\left(n+\frac 12\right)\Omega\right)}{\sin\left(\frac \Omega 2\right)}$ for $(n+1/2) = 10/2$ because $n \in Z$
In DTFT , in my book there is no property like in continous time to transform convolution in $\Omega$ domain to multiplication in time domain so I don't know what to here as well.

update:
After some comments and help from people who answered :
I am going to try to do it as juch $\frac{sin(10\Omega/2)}{sin(\Omega/2)}= \frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}e^{j\Omega(10-1)/2}=\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}\Big]e^{j9\Omega/2}$
I am to take advantage of the property : $\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}\Big]e^{j9\Omega/2} \rightarrow 2\pi F^{-1}{\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}\Big]} * F^{-1}[e^{j9\Omega/2}]$
The result is :
$F^{-1}[e^{j9\Omega/2}] =$
$\frac{1}{2\pi}int_{\pi}^{\pi}e^{j9\Omega/2}e^{j\Omega n}d\Omega = \frac{1}{2\pi}\frac{e^{j\Omega(9/2 +n)}}{j(9/2+n)}\Big|_{-\pi}^{\pi}=\frac{4(-1)^n}{2\pi(n+9)}$ ( i think)
and $F^{-1}[e^{j9\Omega/2}]=1$ for $n \in [0,9]$ and 0 anywhere else.
Now we need to compute the convolution of those 2:
the result should be non - zero only when $n \in [0,9]$ so:
$F^{-1}\Big[\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}\Big]e^{-j\Omega(10-1)/2}\Big] = \begin{cases} \frac{4(-1)^n}{(n+9)} & n \in [0,9] \\
                         0 & else \end{cases}$


texts:
I understand that truncating a signal in time 'smears' the frequency response depending on the window chosen.  In general, the shorter the signal duration, the more 'flattened' the frequency response, as seen here (http://www.thefouriertransform.com/pairs/box.php):

However, how does the window length affect the frequency response of the (bandlimited additive white Gaussian) noise?  Assume a rectangular window of amplitude $A$, duration $T$, and a corresponding $\operatorname{sinc}(\cdot)$ main lobe in the frequency domain with amplitude $A\,T$ and width $\frac{2}{T}$:

$$\begin{align}
\mathscr{F}\bigg\{A \cdot \operatorname{rect}\left(\tfrac{t}{T}\right) \bigg\}
 &= \int_{-\infty}^{+\infty} A \cdot \operatorname{rect}\left(\tfrac{t}{T}\right) \, e^{-j2\pi ft} \ dt \\
\\
 &= \int_{\frac{-T}{2}}^{\frac{+T}{2}} A \, e^{-j2\pi ft} \ dt \\
\\
 &= A \, \frac{\sin(\pi fT)}{\pi f} \\
\\
 &= A \,T \, \operatorname{sinc}(fT) \\
\end{align}$$

If $A$ was fixed, and $T$ was halved, that would result in a $\operatorname{sinc}$ of halved amplitude but doubled main lobe width.  It would then seem that convolving this $\operatorname{sinc}$ would result in the 'same' amplitude of noise in the frequency domain because of the $\frac{1}{2}\cdot 2 = 1$ canceling.  That is, the effective noise bandwidth contributing to a given frequency is doubled, but the contribution per Hz of that bandwidth is halved.

Is that true? And in general, how does a window's duration and shape affect the frequency response of noise?
If (1) is true, does this imply that halving the window duration will also halve the SNR of a single sinusoid?  (Because the sinc magnitude of the signal is halved, but the noise floor remains constant)

Edit: One point I realized is that there may be destructive interference among noise components of different frequencies, and therefore this is not such a simple analysis as just convolving the fourier transform of the window function with the square root of the noise power spectral density.  Perhaps uniformly distributed noise phase at each frequency could be assumed?

I don't have access, but perhaps this paper is useful? http://ieeexplore.ieee.org/document/199437/


texts: Solving a second order difference equation

texts:
I need to solve the following 2nd order difference equation

$$3y[n]-y[n-1]+y[n-2]={\bigg(\frac{1}{2}\bigg)}^{n} u[n]$$

where $u[n]$ is Heaviside's step function. I started by trying to find the homogeneous solution:

$$3-z^{-1}+z^{-2}=3-\frac{1}{z}+\frac{1}{z^2}=3z^2-z+1$$

Homogeneous equation:

$$3z^2-z+1=0\leftrightarrow z=\frac{1\pm\sqrt{1-4\times 3\times 1}}{2\times 3}\leftrightarrow z_1=\frac{1+\sqrt{11}i}{6};z_2=\frac{1-\sqrt{11}i}{6}$$

Since I have complex numbers I don't know how to write the homogeneous solution with the constants $C_1$ and $C_2$. I also don't know how to find the particular solution.


texts: Frequency Axis of Discrete Fourier Transform (DFT) with Odd Number of Data Points

texts:
I am trying to understand the logic behind making a frequency axis in DFT. I am using for time based light absorbance. When we have even number of data points (N= even integer), collected over a length of time A, the lowest frequency $\omega_1$ that can be resolved is 1/A, the next step is 2/A, the $k^{th}$ step is k/A, where k =1,2,3..., N/2. The highest frequency $\omega_{max}$ = N/2A based on Shannon's or Nyquist criterion. This the maximum value on the frequency axis.

How do we scale the frequency axis of FFT when we have odd number of data points (N+1) collected over a length of time A, where N is a even number? I cannot find similar reasoning for scaling when we have odd number of data points except in the book excerpt below.

The step size of frequency in case of odd number of data points:  k/A whether the data points are even or odd.

Please see the excerpt from a book called DFT: An Ownwer's Manual by Briggs. What is the author trying to say when he says that the highest frequency $$\frac{N}{(2A)}$$ does not quite coincide with the endpoint of the frequency domain which have the values $\Omega$/2?Please note that the author is using N as an even number.

Is he indicating that there is a slight error when we have a DFT of an odd number of datapoints?


texts: What do you mean by radio sensing a channel?

texts:
I am slowly getting into the field of wireless communications, in particular with WiFi based communications.

802.11 MAC mentions that

  A station wishing to transmit first listens to the medium for a predetermined period and If it senses to be IDLE , then it transmits . If the medium is BUSY , it waits for some random time  and then again senses the channel

I have two questions here:

What exactly do you mean by SENSING here ?
I have a radio A which wants to transmit on 2.45 GHZ. Do I transmit anything dummy data to know if its is buysy or idle ? I am thinking of programming wise
how does radio A knows that channel is busy or idle ? Assuming there is no ACK concept here ?


texts: How to harvest the highest amount of information from an image?

texts:
I wrote a script that detects objects, crops the part of the image where an object is located and tries to increase the quality of cropped image. The object detection runs tensorflow. To increase image quality I use the Google RAISR AI. It tries to increase the image's resolution by assuming "in-between-pixels" through its algorithm. Although this method is not perfect, it fits my need. After that, a threshold of the image is built for further analyzing.

Original image with detected object (upper line in object = reflection):

Cropped object image:

Improved image with RAISR:

Threshold:

With these images, my script can create a threshold. The main aspect of success is the quality of the object as well as the quality of its reflection on the inner electrode of the pipe.

Now I have datasets which show a much weaker reflection of the object, for example when the object is not at the inner electrode but at the outer wall.

Original image with detected object:

Cropped object image:

Improved image with RAISR:

As you can see, no reflection is visible on this image with the human eye.

This makes it difficult to build a threshold:

The reflections can only be seen on a few images of the dataset. These could be used to build a threshold. Unfortunately 98% of the images are not of this "quality".

I am looking for a method to "harvest" as much information as possible from images like these:

and build a threshold. Are there methods to make the script "see" the reflection of the object although it can't be seen with the human eye? This image should result in this threshold:

I know that my wording of the problem might be a bit unprofessional. I hope I could make it clear.


texts: Does the output of this system have a closed form expression?

texts:
\begin{align}
x(n)&=\left(-1\right)^n u(n)\\
h(n)&= u(n-1)+u(n-2)\\
y(n)&= x(n)*h(n)\\
y(n)&= \sum_{k=-\infty}^\infty \left(-1\right)^k u(k)\left[u(n-k-1)+u(n-k-2)\right]
\end{align}
Graphically it looks like the output is $u(n-1)$ but I can't do the math to get it into the closed form.

Edit: Ok I pretty much worked it through I'm just iffy on justifying $n \ge 1$

$u(k)$ sets the bottom limit of the summation to $0$ and $u(n-k-1)$ sets the upper limit to $n-1$. $n$ has to be at least $1$ for the summation to make sense and the partial sums = $1$.


texts: Trying to create an Extended Kalman Filter for this problem at hand

texts:
Currently I have a system that measures the GPS coordinates of an object. The object is first detected and then using trigonometry, the GPS coordinates are determined, as we know of the GPS coordinates of the camera itself.

However, the camera is on a moving object and therefore the data for the GPS coordinates can be quite noisy. In order to tackle this, I have decided to use an EKF in accordance with the system at hand.

In order to integrate this system into the program that I am currently using, I have decided to use this library here: https://github.com/simondlevy/TinyEKF. However, here it says that

  TinyEKF requires you to write only a single model function, filling in
  the values of the state-transition function f(x), its Jacobian matrix
  F(x), the sensor function h(x), and its Jacobian H(x). The prediction
  and update then handled automatically by passing the observation
  vector z to the step function.

The output of my GPS code, is always [lat, lon]. What is x and what is z here? I don't understand how to create a state transition matrix. however, I am aware that a jacobian can be calculated from the state transition matrix, if I know what it is. Also, what is the sensor function(x)?


texts:
I heard that an IIR filter has delay and phase shift effect, but not clearly what are they. For example, 1-pole Lowpass has 90 degree phase shift, 2-pole Lowpass has 180 degree? So does it mean two 2-pole Lowpass filters has 360 degree phase shift? Beside the HighPass, LowPass, AllPass filter, other IIR filter like LowShelf, HighShelf, Peak... has phase shift too?
When I apply an IIR filter using biquad transfer function, does it has any phase shift diffirent if I use direct form I (4 delay registers) vs direct form II (2 delay registers)?
If IIR filter has phase shift, the phase shift occurred at the filter's target frequency, or all the frequencies? How can I phase shift (all frequencies) 90 degree clockwise & counterclockwise just using AllPass filter?
I know some of my questions are noob, but I'm just study dsp myself, and there're somethings that I really confused about. Sorry for my bad English too :)


texts: QPSK demodulator on GNURadio continuously loses phase lock

texts:
I am designing a QPSK demodulator to be used for receiving Meteor M2 (LEO satellite) signals (QPSK at 72Kbaud). However, I am finding a lot of problems in my way.
First of all, this is the flowchart I came with:

First, I filtered the signal, applied it a $\frac{3}{16}$ fractional downsampling and fed it to a Polyphase Clock Sync block (8sps input, 2sps output and 2pi/100 bandwidth, 0.5 TED Gain and 0.7 Damping factor and 0.15 Deviation). Then, since I thought that for lower passes multipath would be a thing, I added a CMA equalizer after the symbol sync block. Lastly, I placed a Costas Loop just before the soft decoder.
When the elevation of the satellite is above ~50°, I obtain a pretty clear constellation with no decoding errors. Nevertheless, the interesting challenge of designing such demodulator is for it to work at lower elevations. It is very important that it does not lose lock at any moment of the pass, since even the shortest lock loss will produce a black line on the final image.
I am having a lot of trouble with the signal lock at low elevations. I assume that this is caused by a "dirty" signal at the input of the Polyphase Clock Sync block.
• Does my flowchart look good? This is the first digital receiver I design in GNURadio and have very little experience in this field.
• Is there any way (other than trial and error) to optimise Symbol Sync block parameters?
• What technique (filtering, equalization, etc.) may I use in order to solve the locking problems?
Thank you in advance!
Edit: the input signal to the low pass filter has already no doppler shift.
Edit 2: can a PLL Carrier Tracking block help in this situation? The signal seems to be frequency stable, but I do not know whether such block can help avoiding phase locks at the symbol sync stage. Can placing the Costas Loop before symbol sync help?
Edit 3: after some refinement, this is the look of the receiver

Thank you very much @gotchi85 por the help


texts: Given f(x) and its FFT F(u), I need to prove that df/dx = F(u)*2iu(pi)/n

texts:
Given $f(x)$ and its FFT $F(u)$, I need to prove that $df/dx = F(u) * 2iu\pi /n.$

$df/dx = f'(x)$ and $n$ is the number of pixels of the one dimensional image $f$.

I tried to use the convolution theorem using conv($f[1 -1]$) which is the same as $df/dx$ and in the frequency domain it's $F(u) * FFT([1 -1])$ however there sizes don't match.


texts: (For c-FastICA) On covariance and pseudocovariance matrix of a complex random vector

texts:
I am currently studying complex FastICA and the paper says that

  Suppose $\mathbf{s}$ is a $n\times1$ complex random vector. If $\mathbf{s}$ has zero mean, unit variance, and uncorrelated real and imaginary part of equal variances, then $E[\mathbf{s}\mathbf{s}^H]=\mathbf{I}_n$ and $E[\mathbf{s}\mathbf{s}^T]=\mathbf{0}_n$.

I don't quite get how $E[\mathbf{s}\mathbf{s}^H]=\mathbf{I}_n$ and $E[\mathbf{s}\mathbf{s}^T]=\mathbf{0}_n$ come about from the conditions.

We have the covariance matrix as
\begin{align}
\operatorname{cov}(\mathbf{s})
&= E[\mathbf{s}\mathbf{s}^H]-E[\mathbf{s}]E[\mathbf{s}^H] \\
&= E[\mathbf{s}\mathbf{s}^H]-\mathbf{0}_{n\times1}\mathbf{0}_{1\times n}\\
&= E[\mathbf{s}\mathbf{s}^H]\\
\end{align}
and the pseudocovariance
\begin{align}
\operatorname{pcov}(\mathbf{s})
&= E[\mathbf{s}\mathbf{s}^T]-E[\mathbf{s}]E[\mathbf{s}^T] \\
&= E[\mathbf{s}\mathbf{s}^T]-\mathbf{0}_{n\times1}\mathbf{0}_{1\times n}\\
&= E[\mathbf{s}\mathbf{s}^T]\\
\end{align}
I don't quite get how to equate the last line of covariance matrix to identity and the pseucovariance to zero.

If I were to write out the matrix,
\begin{align}
E[\mathbf{s}\mathbf{s}^H]
&=E\left\{\begin{bmatrix}
s_1s_1^* & s_1s_2^* &\cdots & s_1s_n^*\\
s_2s_1^* & s_2s_2^* &\cdots & s_2s_n^*\\
\vdots & \vdots &\ddots & \vdots\\
s_ns_1^* & s_ns_2^* &\cdots & s_ns_n^*\\
\end{bmatrix}\right\}
\end{align}
and
\begin{align}
E[\mathbf{s}\mathbf{s}^T]
&=E\left\{\begin{bmatrix}
s_1s_1 & s_1s_2 &\cdots & s_1s_n\\
s_2s_1 & s_2s_2 &\cdots & s_2s_n\\
\vdots & \vdots &\ddots & \vdots\\
s_ns_1 & s_ns_2 &\cdots & s_ns_n\\
\end{bmatrix}\right\}
\end{align}
I still can't quite figure how all of these eventually becomes identity and zeros.


texts:
I have an MRI K-Space data 320 x 320 x 256 x 8 (4D complex double) from < mridata.org > . The data represents 320 x 320 K-Space of 256 slices from 8 channels. I am trying to reconstruct images for each slice. Here is the matlab code I tried:

% Version - 1
kspacedata= kspacefile.kdata;
imRef = ifftshift(ifftn(kspacedata));
imSOS    = squeeze(sqrt(sum( abs(imRef).^2, 4))); % sum-of-squares to combine all channels

% Version - 2
kspaceSlicedata = kspacedata(:,:,100,:);
imSliceRef = ifftshift(ifftn(kspaceSlicedata));
imSliceSOS    = squeeze(sqrt(sum(abs(imSliceRef).^2, 4))); % sum-of-squares to combine all channels

% Plotting
figure;
subplot(1,2,1); imagesc(imSOS(:,:,100));title('Version - 1');axis image; colormap(gray);
subplot(1,2,2); imagesc(imSliceSOS);title('Version - 2');axis image; colormap(gray);

In Version - 1, I take inverse transform on entire kspace data and plot image data corresponding to slice 100. In Version - 2, I take kspace data corresponding to slice 100 and do inverse transform on this kspace and plot the image data. Here is the output images I obtained.

I thought both versions will reconstruct the same way. But images appeared different. How can I take kspace data corresponding to one slice and perform inverse transform to reconstruct slice image?

Thanks.


texts: How to convert a filtered time-domain signal into frequency one?

texts:
I have a signal in matlab like this:

I have designed a band-pass filter in 16-32 MHz range:

and applied the filter to the signal,

the result is like this :

My Question:

what is the easiest way to convert this time-domain signal into a frequency based signal?

can I export the signal into workspace and use the data variable?


texts:
Let's assume we have an $x(n)$ time sequence, whose $f_s$ sample rate is 20 kHz. We are required to design a linear-phase lowpass FIR filter that will attenuate the undesired high-frequency noise beyond 4kHz analog frequency. So we design a lowpass FIR filter and come out with an equation for the unit impulse response $h_{low}(n)$,and assume that our filter design exercise is complete. Sometime later, unfortunately, we learn that the original $x(n)$ sequence's sample rate was not 20kHz, but 40 kHz. What must we do to our lowpass filter's $h_{low}(n)$ coefficients, originally designed based on a 20kHz sample rate, so that they will attenuate $x(n)$'s undesired high-frequency noise when the $f_s$ sample rate is actually 40kHz?

Hint-typical low pass filter's frequency response makes a pair with its unit impulse response which is a sinc function


texts: Low samples numbers machine vision training for object detection

texts:
I'm trying to set-up a machine vision application in which I have to detect the position of a printing nozzle (in the red sqaure in the image) in a field of view in order to prcisely position a ganty stage.

I currently struggle to detect the wished features. I have tried to implement a Haar cascade using opencv but the detection is not performing.

To train my classifier I used the full image with the location of the nozzles blacked out as negative and the cutout of the red rectangle as positives.
I have few positives and negatives (15 of each) that are basically the same image with different illuminations/expositions and placements of the printer head in the field of view.

Since this approach did not work, what kind of image recognition should I use in that case? Would a training of a custom set in tensorFlow give better results? I've seen implementations of DNN using existing vast libraries of tagged images but I'm not sure how if this is the best method for this application.

Thanks in advance


texts: Compressed Sensing Mathematical Concept in Signal Processing

texts:
I am new in the field of compressive sensing, I've read many papers explaining that compressive sensing is used widely in sparse signal reconstruction. I've tried to understand how compressive sensing is used in signal processing, but I couldn't understand the concept.

According to the paper On some common compressive sensing recovery algorithms and applications - Review paper in page 2, THE MATHEMATICAL BACKGROUND OF THE COMPRESSIVE SENSING CONCEPT, I got the best description for compressing sensing mathematical concept. Assuming we are dealing with signal $x$ with length $N$, that is sparse in time domain. So according to equation 5, we can write the system equation as $Y = AX$, where $A$ is the CS matrix. My questions,

1- What is the signal $x$ with length $N$ which we want to reconstruct ? Is it $X$ or $Y$ ? And what does $A$ represent?

2- Is there any tutorial or file which explain that point with details and in easy way?


texts: Counting the number of groups of 1s in a boolean map of numpy.array

texts:
I am right now dealing with some image processing in Python via PIL (Python Image Library). My main aim is counting the number of colored cells in an immunohistochemistry image. I know that there are relevant programs, libraries, functions and tutorials about it, and I checked almost all of them. My major aim is writing the code manually from scratch, as much as possible. Hence I am trying to avoid using lots of exterior libraries and functions. I have written the most of the program. So here is what's going on step by step:

Program takes in the image file:

And processes it for the red cells (basically, it turns off the RGB values below a certain threshold for red):

And creates the boolean map of it, (gonna paste a part of it since it is big) which  basically just puts 1 wherever it encounters with a red pixel in the processed second image above.

22222222222222222222222222222222222222222
20000000111111110000000000000000000000002
20000000111111110000000000000000000000002
20000000111111110000000000000000000000002
20000000011111100000000000000000001100002
20000000001111100000000000000000011111002
20000000000110000000000000000000011111002
20000000000000000000000000000000111111002
20000000000000000000000000000000111111102
20000000000000000000000000000001111111102
20000000000000000000000000000001111111102
20000000000000000000000000000000111111002
20000000000000000000000000000000010000002
20000000000000000000000000000000000000002
22222222222222222222222222222222222222222

I intentionally generated that frame like thing on the borders with 2s to help me with counting the number of groups of 1s in that boolean map.

My question to you guys is, how come I can efficiently count the number of cells (groups of 1s) in that kind of boolean map? I have found http://en.wikipedia.org/wiki/Connected-component_labeling which look extremely related and similar, but as far as I see, it is at the pixel level. Mine is at the boolean level. Just 1s and 0s.

Thanks a lot.


texts: Interpolating on the borders of differently-resolved images

texts:
I'm creating a three-dimensional model of the earth based on SRTM height data. The data set is pretty huge, so only a small fraction of the data is available at any given time.

The height data is organized in 1° x 1° tiles which come in different resolutions (from 1x1 to 1200x1200, roughly in power-of-two steps). For visualization purposes, I'm sampling the surface in equidistant steps. The sampling theorem gives me an optimal tile resolution to work on, and as long as the point is in the middle of a tile or all neighbouring tiles have the same resolution, this is pretty straightforward.

I'm currently employing cubic splines to get a continuously differentiable surface function. Data points are grey, the sampled point is the red x, used data points are green.

The problem appears when a point is near the border of two differently-resolved tiles.

Using the point's nearest neighbours and linearly interpolating between them gives me a continuous surface, but not a differentiable one:

The most common of such neighbourships is unfortuanely a large (1200x1200 or 600x600) tile being next to a 1x1 one.

I need an algorithm for interpolating smoothly on these borders that offers

near-exact results if the queried point coincides with a tile data point
continuous differentiability in both directions
high-frequency features fading out quickly into low-resolution tiles
good performance (possibly at the expense of accuracy)


texts: How to perform an FFT on a signal with a sampling rate of 44100?

texts:
This is probably a very simple Signal Processing question but I can't quite work it out. I'm currently writing an android app which needs to perform an fft on audio data taken from the phone's microphone. I'm currently reading in the audio data into a byte array buffer at a sampling frequency of 44100Hz.

 int bufferSize = AudioRecord.getMinBufferSize(44100,
                        AudioFormat.CHANNEL_IN_STEREO, AudioFormat.ENCODING_PCM_16BIT);

                final AudioRecord recorder = new AudioRecord(MediaRecorder.AudioSource.MIC, 44100,
                        AudioFormat.CHANNEL_IN_STEREO, AudioFormat.ENCODING_PCM_16BIT,
                        bufferSize);

                final AudioTrack audioTrack = new AudioTrack(
                        AudioManager.STREAM_MUSIC, 44100,
                        AudioFormat.CHANNEL_IN_STEREO, AudioFormat.ENCODING_PCM_16BIT,
                        bufferSize, AudioTrack.MODE_STREAM);

                final byte[] buffer = new byte[bufferSize];
                final int finalBufferSize = bufferSize;

Once the buffer has been filled I convert the data to an array of Complex numbers and pass the data into an fft function. My problem is that the fft function I have cannot work with 44100 samples as it requires the length of the array of samples to be a power of 2. I've found what I think are some solutions online but I'm a little unsure if I'm fully understanding them correctly.

As I understand it I can either pad the array out with zeroes in order to make the array length a power of 2 or I can take a fixed fft size, e.g. 1024 samples. If I use the second solution do I need to divide the original array of samples into 1024 chunks and run an fft on each? Any advice on how I should go about this would be appreciated.

Here is the fft class I'm using:
http://introcs.cs.princeton.edu/java/97data/FFT.java.html


texts: Is It Possible to Reconstruct the Original Image from Laplacian Pyramid?

texts:
In our computer vision course, my professor said we can reconstruct the original image from the Laplacian pyramid. It just feels to me that information is lost after doing subsampling and upsampling.

Is it possible? If so, why? Are the reconstructed image exactly the same as the original pixel by pixel? And is this answer correct?

Generate Laplacian pyramid $(f_2, h_1, h_0)$ from image $f_0$:

Reconstruct $f_0$ from Laplacian pyramid:

source


texts: can $\frac{1}{H(z)}$ be causal and stable?

texts:
if we have linear phase FIR filter $H(z)$ which is causal and stable
can $\frac{1}{H(z)}$ be causal and stable ?

can it be causal without been stable ?


texts: number of possible component in sinusoidal model

texts:
Suppose that we have the following model

$y(t) = A_1\sin(\omega_1 t+\phi_1) + A_2\sin(\omega_2 t+\phi_2) + ... + A_p\sin(\omega_p t+\phi_p) + z(t)$

My question is not related to how to determine the number of components, but I am interested in the maximum. For example, let us assume that we have sampled this model and get sample with size let's say $294$, the question arises what can be the maximum number of components? For instance, can it be more then $147$? or more then half the size of the signal sample? Or can it be in general something $20-50$, for instance? I am not expecting an exact number, but a possible maximum number which can be in this case.


texts: Determine g(t) for minimal symbol error probability, given h(t), K=1 and Minimum distance receiver

texts:
Determine $g(t)$ for minimal symbol error probability, given $h(t)$, $K=1$ and Minimum distance receiver.
$h(t)=   \begin{cases}
      0, & \text{$t<0$}\\
      2, & \text{$t\in(0\;, 10^{-3})$}\\
      1, & \text{$t\in(10^{-3}\;,2\cdot 10^{-3})$}\\
      0, & \text{$t>2\cdot 10^{-3}$}
    \end{cases}  $
This is ISI-free if $h(t)*g(t)|_{t=kT}=
    \begin{cases}
      1, & \text{if $k=0$}\\
      0, & \text{if $k\neq0$}
    \end{cases}  $
Using a matched filter as $g(t)=ch^*(-t)$ the convolution becomes
$$(g*h)(t)=\int_{-\infty}^\infty{h(\tau)g(-\tau)d\tau }=c\int_{-\infty}^\infty h(\tau)h^*(\tau)d\tau=1$$
How do I calculate the $\int_{-\infty}^\infty h(\tau)h^*(\tau)d\tau=1$ integral or the pulse energy?
From what I understand it is the same as the pulse energy (which I don't know how to calculate either).
From here I understand how to determine $g(t)$ and it will be the $c$-scaled and time inversed variant of $h(t)$, but I don't really understand how to calculate the pulse energy.
Is $E_h=\int_{-\infty}^{\infty}|h(\tau)|^2d\tau$$\rightarrow$
$E_h=\int_{0}^{10^{-3}}2^2d\tau+\int_{10^{-3}}^{2*10^{-3}}1^2d\tau$ the correct way?


texts: Signal Reconstruction in Compressed Sensing with a Simple Vector Signal as an Example

texts:
While going through the different types of reconstruction algorithm as mentioned in Richard G. Baraniuk - Compressive Sensing - Lecture Notes (Also on DocDroid), I came to know that minimum $ {L}_{1} $ norm reconstruction is better than minimum $ {L}_{0} $ and $ {L}_{2} $ norm reconstruction.

Can anyone show me why $ {L}_{1} $ norm is better by taking a simple signal vector of small length as an example?
I would like to know what will be the final optimization results in each reconstruction method.


texts: Formula of frequency modulated sine wave from data points

texts:
I'm trying to determine the formula of a variable frequency sine wave based on the values of max1, max2, min1 and min2. y(0) isn't necessarily 0 and I can assume the frequency increases or decreases linearly (I only want to plot the first couple oscillations).

I've figured the formula for the non-modulated signal (based only on max1 and min1) as the following:

$y(x) = sin(2π*f*x + C)$
where
$f=\frac{1}{2\left(\max _1-\min _1\right)}$
and
$C = -\left(max_1-\frac{\left(max_1-min_1\right)}{2}\right)\cdot \frac{\pi }{\left(max_1-min_1\right)}$

I've also figured that plotting:

$y(x) = sin(2π*(f+M(x))*x)$ where
$M(x)= a*x + b$

results in an increasing frequency signal but I can't figure out the exact formula to have the wave match the data points.

Any help would be much appreciated, my math skills are rusty...


texts: Approximating Lorentzian Fourier transform with FFT

texts:
I have a Lorentzian frequency distribution

$F(w) = \frac{1+iz}{1+z^2}$

Where

$z = \frac{w-\Omega}{R}$

With $\Omega$ being the peak frequency and R the decay constant. I know that analytically the Fourier transform should be

$F(t) = exp(i\Omega 2\pi t)exp(-Rt)$

When i take the FFT of this expression, it doesn't return the original frequency plot. I understand that there may be a scaling factor (1/n) that should be in there somewhere, but even when i scale for one frequency, if i then change $\Omega$ or R, the amplitude is no longer scaled properly, suggesting that the scaling factor is a function of $\Omega$ and/or R. The FFT also seems to be mirrored along the intensity axis.

I'm fairly new to DSP, but i do understand that the continuous fourier transform is not the discrete fourier transform. I've read this (https://dspillustrations.com/pages/posts/misc/approximating-the-fourier-transform-with-dft.html) but that approach makes the approximation worse.

I'd like to return the original frequency distribution when i take the FFT of my time signal. Am i missing something fundamental or is it a fairly simple scaling error? I've attached my code below.

Cheers.

# R script to compare FFT and Analytical fourier transform
library(SynchWave)

#-------------------------------------------------
# Frequency and time axes
n <- 100
f <- seq(0, 1, length.out = n)
t <- seq(0, n, length.out = n)

# peak paramaters
O <- 0.3 # Frequency values from 0->1
R <- 0.04 # Decay in arbritrary units

z <- (f-O)/R

# The original lorentzian frequency
ff <-complex(re = 1, im = z)/(1 + z^2)

# creating the time domain signal
ftideal <- exp(-R*t)*exp(complex(i = (O)*2*pi*t))

unscaled <- (fft(ftideal))
scaled <- unscaled - min(Re(unscaled))

plot(f, Re(ff), type = 'l')
lines(f, Re(scaled), type = "l", col = 'red')

```


texts: Hermitian FFT vs discrete Fourier transform. Python's $\tt hfft$ vs $\tt fft$

texts:
What is the difference between the Hermitian FFT and DFT? Particularly in Python, there are two functions fft and hfft.

numpy.fft.hfft(signal) vs numpy.fft.fft(signal)

What I simply could find out is:
The Hermitian has to do something with symmetry and needs 50 times longer to calculate, while producing a 'slightly' different result than the 'discrete' FFT. (tested on an audio file of machinery sounds and length of 1.5 sec).


texts: How to Save FFT output for fast read-in?

texts: For a project, I'm reading in hundreds of short (4 sec) wavefiles and giving them to scipy fft for the frequency information for further processing/experimenting in various ways. Problem is, I want to run the program many times while changing params, code, test new ideas, etc, but each time I'm rereading, re-fft-ing the original wavefiles. I'm wondering if there is a smart way to "save" the fft numpy arrays for quick retrieval, so I'd only have to read them in once, if this is possible and not a dumb idea.


texts:
What techniques/algorithms can I use to remove noise from a raw recording of sound (voice)?

The purpose is to get a smoother graph (removing the "jaggedness"). What I have tried was to average-out small deviations by using the surrounding two values to either side, but ended up with distorted sound. What better way is there to "smooth-out" the graph?


texts: The optimal estimator using FIR Wiener filter

texts:
The received signal $y(n)$ is defined as
$$y(n) = h(n)*x(n) +v(n)$$
where $h(n)$ is unknown impulse response and $v(n)$ is a white Gaussian noise.
The value of $x(n)$ for $n=1, ..., 32$ are given. The task is to calculate an estimator of $x(n)$. I calculated the estimator using FIR Wiener filter which gives the optimal means error, when the impulse response $h_{wie} = R_y^{-1} r_{xy}$ where $R_y=E[ YY^T]$ and $r_{xy}= E[xy]$. My intuition is when we have more value of $x(n)$, we should get the better estimator but that not what i can observe from my computation i get optimal estimator when filter orders is around 24-30 by try and error.
Can someone explain a systematic way to find optimal taps/orders of FIR Wiener filter? And is there a way to improve the estimator i.e different types of filters?


texts:
How does (if it does) frequency affect/relates to Data Rate.

Why for example using the same modulation, LoRa (in the physical layer),
you get somewhere around 250kbps at 2.4Ghz and only 20kbps at 868Hz ?


texts: Averaging frequency spectra gives incorrect RMS

texts:
I am averaging several FFT results together in an attempt to get the average frequency spectrum of a track of audio. I'm attempting to use Parseval's Relation to get the RMS of the audio from the resulting average spectrum, but i've noticed the RMS i get is wrong, usually only by a couple dB. Why does this happen? Does it mean i am averaging the spectra together incorrectly? Is it more complicated than simply adding the corresponding N magnitudes together and dividing by N? When i get the RMS from a single FFT using Parseval's Relation it is correct.

Here is the code i am using to average the FFTs:

    int N = pow(2, resolution);

    //temporary buffers
    vector< array<double,2> > RXT;
    vector< array<double,2> > IXT;

    //final signals
    REX.assign(N, array<double,2>() );
    IMX.assign(N, array<double,2>() );
    magnitude.assign( N, array<double,2>() );
    phase.assign( N, array<double,2>() );

    int start;
    bool lastSeg = false;
    int seg = 0;
    //double radians = 2 * PI / N; //used for windowing functions

    while ( !lastSeg ) {

        start = seg * N;

        if (X.size() >= start + N) {
            RXT.assign(X.begin()+start,X.begin()+start+N);
        } else {
            RXT.assign(X.begin()+start,X.end());
            for ( int n = 0; n < RXT.size(); n++ ) {
                for ( int chan = 0; chan < channels; chan++ ) {

                    RXT.at(n).at(chan) *= (double)N / RXT.size();

                }
            }
            RXT.resize( N, array<double,2>() );
            lastSeg = true;
        }

        IXT.assign( N, array<double,2>() );

        FFT(RXT, IXT);

        for ( int n = 0; n < N; n++ ) {
            for ( int chan = 0; chan < channels; chan++ ) {

                magnitude.at(n).at(chan) += getMag( RXT.at(n).at(chan), IXT.at(n).at(chan) );

            }
        }

        seg++;

    }

    for ( int n = 0; n < N; n++ ) {
        for ( int chan = 0; chan < channels; chan++ ) {

            magnitude.at(n).at(chan) /= seg;

            REX.at(n).at(chan) = getReal(magnitude.at(n).at(chan), phase.at(n).at(chan));
            IMX.at(n).at(chan) = getImag(magnitude.at(n).at(chan), phase.at(n).at(chan));
        }
    }


texts: confusion regarding correlation in noise?

texts:
I am reading online article about image processing at following link

https://www.electronicsforu.com/electronics-projects/software-projects-ideas/image-processing-using-matlab-basic-operations-part-2-4

I have also attached a snapshot of a paragraph of this link and i have highlighted confusing statements

As paragraph says

"In either case, the noise at different pixels can either be correlated or uncorrelated".

That give rise to two queries

1) What is difference between correlated noise and non correlated noise?

2)If an image has noise, all pixels will have same noise,correlated or non correlated, or is it possible that some pixels have correlated noise and others have non correlated noise?


texts: Filter ECG signal with Wavelet and Python

texts:
I have to filter the signal of an ECG with the wavelet method with Python. The data is in a txt file. I can create my dataframe with pandas, display that with seaborn, but can not find a way to apply the filter.

The wavelet method is imposed. I do not really know how to do it.

First I tried to understand the mathematical formula to transcribe it into a python algorithm: it's clearly beyond my abilities and my knowledge (and filling this gap would take me months, even years)

I also tried using the python pywavelets library.

The problem is that I do not know what to do with cA and cD :

cA, cD = pywt.dwt([1, 2, 3, 4], 'db1')

My research on google (in English and French) directs me towards papers of scientific research which are also outside my capacities of comprehension.


texts: Uplink Pilot Design in TDD based OFDM system

texts:
In a paper “Design of Orthogonal Uplink Pilot Sequences for TDD Massive MIMO under Pilot Contamination”,it was mentioned that, In the uplink channel estimation, if each user occupy all available subcarriers, it will not be suitable for the TDD massive MIMO systems. Because in TDD massive MIMO system only a limited portion of the coherence interval can be used for the uplink training.
In this, I don’t really understand why coherence interval comes into picture while taking about subcarriers. Also, why it is not suitable if every users use all available subcarriers?
Thanks.


texts: How do I zero pad time samples if I only have frequency samples?

texts:
Related to this question:
Convolution in frequency space with fft and ifft -> num output samples?

When doing convolution in time domain by doing multiplication in frequency domain, you apparently can zero pad the time domain samples if you want the time domain convolution result (ifft of the fft multiplication result), to make the resulting time domain convolution be the "correct" length (len(a)+len(b)-1).

What I'm wondering is what do you do if you don't have the time domain samples yet?

I'm generating two sets of IFFT bin data, doing convolution by multiplying in the frequency domain, and then using IFFT to turn the convolution result into time domain data.

When i do this, the result is the same length as my IFFT bins though, instead of the length i would expect.

How do I lengthen the result in the frequency domain?

Thanks!


texts: Confusion understanding digital signal processor?

texts:
I am reading wikipedia article about digital signal processor and as shown in end of attached snapshot
." Most general-purpose microprocessors can also execute digital signal processing algorithms successfully, but may not be able to keep up with such processing continuously in real-time"
Please elaborate with example/application? What is meant by real-time here in context of dsp?


texts:
My system is described by the following ODE:

$$
\frac{\mathrm{d}y(t)}{\mathrm{d}t} = a-y(t)x(t)
$$

Where $a$ is a constant and x(t) is a poisson process so that:

$$
E[x(t)x(s)] = qI\delta(t-s)
$$

The solution is of the form:

$$
y\left( t \right) = e^{-\int^t{x(z)dz}}{\int^t{a\,e^{\int^z{x(w)dw}}}dz+ c}
$$

I would like to find a close form for the output total power $E[y^2(t)]$ in function of time.


texts: How to recover frequency and damping coefficient from superposition of damped oscillators?

texts:
I have a signal which is (more or less) the superposition of multiple damped oscillators.  There are relatively few oscillators involved, 2-3 would account for most of the signal.  They have periods which are not multiples of each other (ie not harmonics), and they all have different damping coefficients / time constants.  They also are not necessarily in phase initially.

I would like to recover the periods, damping coefficients and phases of the oscillators from the signal.

A synthetic example of a signal like this can be produced by the following python code:

from pylab import *
x = linspace(0,1000,num=1000)
y1 = 1 * sin(2*pi*x/11 + 0) * exp(-x/300)
y2 = 0.5 * sin(2*pi*x/17 + 1.2) * exp(-x/300)
y3 = 2 * sin(2*pi*x/50 - 0.25) * exp(-x/50)
y = y1 + y2 + y3
figure(figsize=(20,10))
plot(y)

This looks pretty similar to the real signal qualitatively, including having some beat frequencies, and a pretty large initial low-frequency component.

This seems like a pretty common problem when analyzing a physical system with multiple coupled oscillators, eg a vibrating object with multiple resonances.

UPDATE Using direct curve-fitting works pretty well, even for strongly damped signals - depending on the choice of initial conditions, that is:

from scipy.optimize import curve_fit

def func(x, a, b, c, d):
     return a * sin(2*pi*x/b - c) * exp(-x/d)

p, pcov = curve_fit(func, x, y, p0=(1,100,0,1000))
print(p)

# [  1.78453351  51.53572084   0.05566558  55.10643076]

This recovers the largest-amplitude (but very strongly damped: time constant = 1 period) part of the signal.


texts: PSS convertion from Frequency to time domain

texts:
Hi I have a question regarding primary sync signal in LTE.

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

-Thanks


texts: How to check if Fourier components are in phase of out of phase?

texts: I have a time series (of 23 timestamps) of which I take the Fourier transform. Now the fourier transform has 23 imaginary values and each has an amplitude and a phase. When I get the phase angle, it given me arctan value of the angle. Now how can I determine if two fourier coefficients are in phase or out of phase.


texts:
I am still relatively new to signal processing and MATLAB and am working on a short-time fourier transform of some data using a prexisting code. Within the code a Hamming window is produced and multiplied by the signal prior to an fft.
The function itself works correctly, however I want the peak of the Hamming window to be over the first sample, rather than the midpoint (i.e shifting the Hamming window back by (no of samples/2)). What would be the simplest way to achieve this zero meaning?
This is the current generation and application of the Hamming window, where ncol is the number of time windows being used in the stft.
y = window(:,ones(1,ncol)).*y;

Any help is greatly appreciated.


texts:
I am reading a nice book on DSP

  The Scientist and Engineer's Guide to Digital Signal Processing

It mentions

  In electronics, the mean is commonly called the DC (direct current) value.

How can a signal have no DC component if the mean is the DC component? Surely every signal's values can be averaged?


texts:
I designed a 1st order Delta Sigma converter in PSpice, except for the decimation filter. I want to implement said decimation filter in Python, rather than in PSpice. Being new to both Python and filters in general, I turned to google, which offered me this "simple" solution.

The explanation is quite clear (I think). In my case, the parameters are the following:
- signal going into the DeltaSigma ADC: 1 kHz sin wave
- sampling rate of the DeltaSigma ADC: 5 MHz
- simulation step size: 100 ns (might be relevant)

When I design my filter (a Windowed Sinc filter, as in the example above) I use the following parameters (I have also used some variations of these parameters):
- cutoff frequency fc = 3 kHz
- transition bandwidth 1 kHz, which results in N = 200001

NB: the filter might be very ambitious, but as far as I can see this is just a computation-time-related problem.

Now, if I implement the above in Python everything seems to work fine up to the filter application to my signal; this means that the sinc filter "looks" nice in the plot, the input signal is ok, etc...

For reference: this is my implementation of the filter.

# Define the parameters
fc = 0.0006
b = 0.0002
N = int(np.ceil((4 / b)))
if not N % 2: N += 1  # Make sure that N is odd.
n = np.arange(N)

# Compute sinc filter.
h = np.sinc(2 * fc * (n - (N - 1) / 2.))

# Compute Blackman window.
w = np.blackman(N)

# Multiply sinc filter with window.
h = h * w

# Normalize to get unity gain.
h_unity = h / np.sum(h)

After this, I simply do a convolution with:

filtered_data = np.convolve(out, h, 'valid')

Image of the output data after the filter:

I would expect the filtered data to be a sin wave, not that "thing".

My guess is that I am missing something fundamental and that somehow I should take care of the fact that the sampled signal is actually sampled with both the 5Mhz clock AND with the 100 ns simulation time.

Any help will be appreciated!

EDIT: clarification on the out signal.

The out signal is generated by a Delta Sigma converter, so supposedly it should be a stream of '1' and '0'.

But since the PSpice simulator has an internal step size (which is variable) the actual values of the out signal have also intermediate values (remember: the output waveform in PSpice is a voltage, so in the transition between 0/1 the voltage is calculated).

So, my out signal looks like this (not actual values, just for clarity):

Time    Out
0       0.000
0.1     0.000
0.145   0.123
0.167   0.679
0.2     1.000
0.3     1.000
0.4     1.000
0.45    0.875
0.493   0.439
0.560   0.293
0.6     0.000

This most likely introduces spurious frequencies in my spectrum.

Working on a possible solution, will let you know if it works!


texts: How to get coordinates(x(time),y(frequency)) of spectral peaks from spectrogram?

texts:
I have created this spectrogram from a wav file. Please have a look:

As it can be observed clearly that my somewhat spectral peaks are visible. I want to get the x = time, y= frequency of these spectral peaks in order to find temporal displacement between the peaks. How can I get these coordinates?

I am looking for a programmatic approach (using python,java etc).

Regards,

Khubaib


texts: If my sampling rate is not high enough to capture all frequencies, can I accurately capture low frequencies?

texts: Let us say that my sampling rate is 1000 Hz. This enables me to accurately capture every frequency in the 0-500 Hz range according to the Sampling Theorem using a DFT. However, there are higher frequencies present in the signal in between 1500-2500 Hz. If I am to convert the 1000 Hz data to its DFT, I understand I will not get the higher frequencies. However, I only care about the lower frequencies. So my question is, will I accurately get the lower frequencies of the signal and their precise amplitudes using the DFT?


texts: How do you draw decision boundaries given this type of 16-qam constellation diagram, I know that there should be no vacant decision boundary but I can't seem to get it.Also, assuming that the correct decision boundary is drawn, how do you derive its total probability of error, if we let the SNR ratios of all symbols are equal?


texts:
Consider the following two statements:

In time axis:

A signal without a DC component is a signal which doesn't have the zero frequency (the DC frequency)
A signal without a DC component is averaged to zero

My questions are:

Why is it incorrect? It sounds like the definition to me.
What is the explanation? It has something to do with the Fourier representation of the signal, I suspect.


texts: Signal processing on hourly atmospheric temperature, wind speed and humidity data for solar panel maintenance

texts:
I am working on a machine learning problem, where i have to predict if a device connected with solar panel at a specific block has been damaged. I have many factors to consider, out of which atmospheric temperature, wind speed and a humidity are to be considered. These 3 factors are being recorded from 3 separate devices. This device helps pack cells ini case of severe windspeed or other environmental factors. What if this device is also damaged. Its to model this.
What sort of signal processing can be done on this data? Do i need to use fourier transformation, wavelet, noise removal or filtering? I am at beginner level in signal processing. I tried doing some research, however i feel that i do need some help in this matter.
Edit: Its not exactly for panel glass if you thinking about the actual panel, but a device that connects between panel and transmission. This device makes decision if panel should be folded when windspeed is high for damage prevention. What if that device has been damaged.


texts: Terminologies - lags, order in time series model

texts:
I am facing some difficulties with the terminologies - lag ($p$) and sequence length (number of data points) (N) used in time series model such as Moving average and Autoregressive model. Considering a Moving Average model described by the equation

$$y(n) = \sum_{i=1}^{p} h(i) x(i-n)$$
 If the number of bits in the input is N=16 defined by $\mathbf{x} = \{1,0,1,1,0,1,0,1,1,0,1,1,0,0,1,1\}$ and the number of lags are p=10
then

Q1: should the model generate a time series of length 'N=16` i.e, would the output of the above model $\mathbf{y} = [y_1,y_2,\ldots,y_N]$ contain 16 elements where $n = 1,2,\ldots,16$?

I am writing out the full expression. I can write the MA process in vector form as $y_n = \mathbf{h}^T \mathbf{x}_n$ where $\mathbf{x}_n = [x_n,x_{n-1},\ldots,x_{n-p+1}]^T$ is the input signal vector and $\mathbf{h} = [h_0,h_1,\ldots,h_{p-1}]^T$ is the MA system coefficients. For example, in programming using Matlab the model would look like for

h = [1,0.2,0.3]; %channel coefficients
N=16;%number of data points
x = rand(1,N); %input data vector
p=2;
y(1) =0.0;
y(2) = 0.0;
for n = 3:N
y(n) = h(1)*x(n) +h(2)*x(n-1)+h(3)*x(n-2);
end

Since the order p=2 would the physical interpretation be that the channel can take in only 2 data values as input? Also, does p=2 mean that the channel has 2 paths through which the data input travels?

Q2: Is there a relationship between the number of input elements and the number of lags?

Please correct me if I have misinterpreted the concepts. Thank you.


texts: Optimizing Two Parallel IIR Filter Sweeps

texts:
Got an audio related question here; So I have two IIR low pass filters that work in parallel to do a filter sweep. Basically:

filtered_sample = filterSweepA(xn) + filterSweepB(xn);

End result: it sounds beautiful! However, now I want to optimize it. I was thinking about combining the two IIR filters into one so I'd only have to call filterSweep only once but I'm pretty sure you can't do that because of the recursive properties...

I think I understand that FIR filters would be the way to go but I've searched the internets and I've seen nothing about FIR filter sweeps.

Is this possible/just difficult to do?
If this can be implemented how would I be doing this?

I don't really know how to figure out FIR filter coefficients so I guess that would be the place to start...

FYI I am controlling the start, stop filter frequency cutoffs and the sweep time (usually around 1.5 seconds).


texts: I know what is angle modulation. I want to know why is named as 'Angle'. Is it because in the vector plot of sidebands are described by various angles? or what?


texts:
I have a same noise passing through two different ARMA filter. The autospectrum of first filter is given by S11=G1*G1*S_in where S_in is input spectrum created by set of random numbers.

Similarly the autospectrum of second filter is given by S22=G2*G2*S_in where S_in is input spectrum created by same set of random numbers as the previous one.

The cross spectrum that i expect is S12=G1*G2*S_in.But the amplitude is reduced is this due to phase shift of the filter or the derivation above is incorrect.

Attached is an image1 explaining the above statements.


texts: How to do edge detection using Curvelet transform?

texts:
I wanted to use Curvelet thresholding for edge detection, that is first take discrete curvelet transform then choose coefficient which are for curve and then reconstruct back to get an edge map.

Please help me on how to do Curvelet thresholding.


texts: IIR Filter - calculating the phase response

texts:
Using formulae from Audio EQ Cookbook from here it was easy to implement biquad IIR filters in C++, i.e. to calculate coefficients b0, b1, b2, a0, a1, a2.

Now, I want to calculate the phase response for a given IIR filter. I have read a number of papers about IIR design, read the websites (DSP Guide) but could not find the way to do it. How I can calculate and specify the phase of an IIR filter during the design?


texts: Refer good books for the following topics in Image processing

texts:
Thanks for reading my question. Please refer good books(links, if possible) for the following topics in Image processing. Your help is very much appreciated.

UNIT II SEGMENTATION

Fuzzy clustering, Watershed algorithm, Active contour methods, Texture feature based segmentation, Model based segmentation, Atlas based segmentation, Wavelet based Segmentation methods.

UNIT III FEATURE EXTRACTION

Texture descriptors- Autocorrelation, Co-occurrence features, Runlength features, Fractal model based features, Gabor filter, wavelet features.

UNIT IV REGISTRATION AND IMAGE FUSION

Image Fusion-Overview of image fusion, pixel fusion, Multiresolution based fusion discrete wavelet transform, Curvelet transform. Region based fusion.